🧪 Historia powstania zbioru danych QM9¶
Zbiór danych QM9 jest jednym z najważniejszych i najczęściej wykorzystywanych benchmarków w chemoinformatyce, uczeniu maszynowym dla chemii oraz projektowaniu cząsteczek. Jego stworzenie było odpowiedzią na potrzebę posiadania dużego, jednorodnego i wysokiej jakości zbioru danych kwantowo-chemicznych, który umożliwiłby rozwój algorytmów ML w obszarze właściwości molekuł.
📜 Geneza projektu¶
Zbiór QM9 został opracowany w latach 2014–2017 przez zespół naukowców z University of Basel, Cornell University i University of Vienna: Ramakrishnan, Dral, Rupp i von Lilienfeld. Został po raz pierwszy opisany w publikacji:
“Quantum chemistry structures and properties of 134 kilo molecules” (Scientific Data, 2014)
Celem autorów było stworzenie systematycznego i kompletnego zbioru obejmującego tysiące małych cząsteczek organicznych oraz ich właściwości obliczonych metodami chemii kwantowej.
🧩 Jak powstał zbiór?¶
Proces tworzenia QM9 obejmował kilka etapów:
1. Wybór cząsteczek (enumeracja strukturalna)¶
Autorzy wybrali wszystkie stabilne cząsteczki organiczne o wzorze sumarycznym zawierającym tylko:
- C,
- H,
- O,
- N,
- F, oraz maksymalnie 9 atomów ciężkich. Dlatego nazwa datasetu to QM9.
📌 W sumie zebrano 133 885 struktur.
2. Optymalizacja geometrii¶
Dla każdej cząsteczki wykonano obliczenia metodą DFT (Density Functional Theory) z funkcjonałem:
- B3LYP
- bazą 6-31G(2df,p)
Była to ogromna kampania obliczeniowa — w sumie miliony iteracji obliczeń kwantowych.
3. Wyznaczenie właściwości fizykochemicznych¶
Dla każdej cząsteczki obliczono ≥ 17 właściwości, m.in.:
- energie atomizacji,
- moment dipolowy,
- HOMO, LUMO, gap,
- częstotliwości drgań,
- polaryzowalność,
- masę,
- geometrię 3D,
- temperaturę wrzenia (przybliżoną metodami kwantowymi).
📊 Dzięki temu zbiór jest idealny do regresji, predykcji właściwości i zadań ML w chemii.
4. Standaryzacja i walidacja¶
Każdą strukturę zweryfikowano:
- poprawność SMILES,
- poprawność minimalnej energii,
- brak fragmentacji,
- poprawna topologia.
Wczesne wersje datasetu zawierały ~300 błędnych geometrii — zostały później usunięte.
🚀 Znaczenie QM9 dla współczesnego ML w chemii¶
Zbiór QM9 stał się złotym standardem w testowaniu modeli:
- GNN (Graph Neural Networks),
- Message Passing Neural Networks,
- SchNet, PhysNet, DimeNet,
- Transformerów molekularnych,
- modeli generatywnych (VAE, diffusion models).
Ponieważ dane są jednorodne, czyste i precyzyjne, QM9 pozwala porównywać algorytmy w sposób uczciwy i powtarzalny.
⭐ Dlaczego QM9 jest tak ważny?¶
- zawiera największy do tamtej pory publiczny zbiór obliczeń DFT,
- jest kompletny w zakresie małych cząsteczek organicznych,
- ma wysoką jakość,
- idealnie nadaje się do projektów ML w chemii ❤️🔥.
Dziś stanowi fundament ogromnej części badań nad AI for Molecules.
📂 Opis zbioru danych QM9¶
W tej sekcji przedstawiono szczegółowy opis struktur danych dostępnych w roboczej wersji zbioru QM9, zawierającej 1000 losowo wybranych cząsteczek. Dane pochodzą z pełnego zbioru 133 885 struktur obliczonych metodami DFT.
🧬 Struktura tabeli danych¶
Po wczytaniu pliku dataset zawiera następujące typowe kolumny (nazwy mogą różnić się zależnie od wersji):
| Kolumna | Opis | Typ |
|---|---|---|
| smiles | SMILES cząsteczki (struktura 2D zakodowana tekstowo) | string |
| formula | Wzór sumaryczny cząsteczki | string |
| rotational_constants_A/B/C | Stałe rotacyjne (GHz) wyliczone z geometrii 3D | float |
| dipole_moment | Moment dipolowy (Debye) | float |
| isotropic_polarizability | Izotropowa polaryzowalność (Bohr³) | float |
| homo | Energia HOMO (eV) | float |
| lumo | Energia LUMO (eV) | float |
| gap | Przerwa energetyczna HOMO–LUMO (eV) | float |
| electronic_spatial_extent | Objętość elektroniczna ⟨r²⟩ (Bohr²) | float |
| zero_point_energy | Energia zerowa (ZPE, Hartree) | float |
| internal_energy_0K | Energia wewnętrzna w 0 K (Hartree) | float |
| internal_energy_298K | Energia w 298.15 K (Hartree) | float |
| enthalpy_298K | Entalpia w 298.15 K (Hartree) | float |
| free_energy_298K | Energia swobodna (Gibbs) w 298.15 K (Hartree) | float |
| heat_capacity | Ciepło właściwe Cv (cal/mol·K) | float |
| mulliken_charges (opcjonalne) | Rozmieszczenie ładunku Mullikena | float/list |
| coordinates_X/Y/Z (opcjonalne) | Współrzędne 3D atomów (Å) | list |
| atom_types | Lista atomów (C/H/O/N/F) | list |
| num_atoms | Liczba atomów | int |
| num_heavy_atoms | Liczba atomów ciężkich (C/O/N/F) | int |
| molecule_id | Unikatowy identyfikator cząsteczki | int |
⚠️ Uwaga: Dokładne kolumny w Twoim Excelu zależą od wersji konwertowanego pliku QM9, ale powyższa lista obejmuje cały standardowy zakres używany w chemoinformatyce.
🧠 Kategorie cech w zbiorze¶
Aby lepiej zrozumieć strukturę danych, można podzielić kolumny QM9 na 4 główne grupy:
1️⃣ Cechy strukturalne¶
- SMILES
- Wzór sumaryczny
- Liczba atomów (ogółem i ciężkich)
- Współrzędne 3D
- Typy atomów
👉 używane głównie do modelowania GNN i obliczeń geometrycznych.
2️⃣ Cechy energetyczne (DFT)¶
- ZPE
- Energie wewnętrzne
- Enthalpia
- Energia swobodna
- HOMO / LUMO / GAP
🔥 kluczowe dla modeli predykcyjnych właściwości kwantowych.
3️⃣ Cechy spektralne i rotacyjne¶
- Stałe rotacyjne A/B/C
- Ciepło właściwe
- Przestrzenne rozciągnięcie ładunku (⟨r²⟩)
4️⃣ Cechy elektrostatyczne¶
- Moment dipolowy
- Polaryzowalność
- Potencjały elektronowe i ładunki (jeśli dostępne)
🎯 Dlaczego te dane są wartościowe?¶
- Pochodzą z precyzyjnych obliczeń kwantowych, nie eksperymentów.
- Charakteryzują się wysoką jednorodnością metodologiczną.
- Idealnie nadają się do trenowania modeli ML predykcji właściwości molekularnych.
- Zawierają różnorodność strukturalną przy ograniczonej wielkości cząsteczek — kluczowe do nauki modeli GNN.
🎯 Cel projektu: Eksploracyjna analiza danych (EDA) zbioru QM9¶
Celem tego projektu jest przeprowadzenie kompleksowej eksploracyjnej analizy danych (EDA) na zbiorze QM9, zawierającym struktury i właściwości kwantowo-chemiczne małych cząsteczek organicznych.
🔹 Główne założenia projektu¶
Zrozumienie struktury danych
- Identyfikacja kolumn, typów danych, braków danych i wyjątków.
- Analiza podziału danych na cechy strukturalne, energetyczne, elektrostatyczne i rotacyjne.
Analiza statystyczna właściwości molekuł
- Rozkłady podstawowych parametrów: energie, moment dipolowy, HOMO/LUMO, GAP.
- Wykrywanie wartości odstających i potencjalnych błędów w danych.
Wizualizacja danych chemicznych
- Graficzne przedstawienie relacji między właściwościami.
- Wizualizacja cząsteczek na podstawie SMILES i geometrii 3D (jeśli środowisko pozwoli).
Tworzenie podzbiorów roboczych
- Wydzielenie mniejszych próbek danych (np. 1000 wierszy) w celu szybkiej pracy i testowania modeli.
Przygotowanie danych do modeli ML
- Identyfikacja cech predykcyjnych i docelowych.
- Wstępne przekształcenia danych i normalizacja dla przyszłych zadań predykcyjnych.
🌟 Wartość projektu¶
- Projekt pozwala pokazać umiejętności w EDA i chemoinformatyce, w tym analizę dużych zbiorów danych chemicznych.
- Prezentuje profesjonalny workflow od wczytania danych, przez ich eksplorację, aż po wizualizację i przygotowanie do ML.
- Stanowi doskonały element portfolio dla kandydatów na stanowiska Data Scientist / Chemoinformatyk / Machine Learning Engineer w chemii.
import pandas as pd
# Ścieżka do dużego pliku
data = r"C:\Users\slast\PYTHON\0_projekty do portfolio\02_EDA_chem\qm9_dataset.xlsx"
# Wczytujemy z opcją ograniczenia do wybranych kolumn lub bez – jak wolisz
df = pd.read_excel(data)
# Losowa próbka 1000 wierszy z zachowaniem miksu danych
sample_df = df.sample(n=1000, random_state=42)
# Zapis nowego pliku roboczego
output_path = r"C:\Users\slast\PYTHON\0_projekty do portfolio\02_EDA_chem\qm9_sample_1000.xlsx"
sample_df.to_excel(output_path, index=False)
print("✔️ Utworzono plik roboczy z 1000 wierszy:")
print(output_path)
✔️ Utworzono plik roboczy z 1000 wierszy: C:\Users\slast\PYTHON\0_projekty do portfolio\02_EDA_chem\qm9_sample_1000.xlsx
duzo recordów tworzenie roboczego df z 1000 wierszami¶
import pandas as pd
data = r"C:\Users\slast\PYTHON\0_projekty do portfolio\02_EDA_chem\qm9_sample_1000.xlsx"
df = pd.read_excel(data)
print(df.head())
X y1 y2 \
0 <rdkit.Chem.rdchem.Mol object at 0x000001BF0ED... -1.291747 1.802935
1 <rdkit.Chem.rdchem.Mol object at 0x000001BEFC9... -0.644026 -0.218808
2 <rdkit.Chem.rdchem.Mol object at 0x000001BF0E0... -0.459969 0.297628
3 <rdkit.Chem.rdchem.Mol object at 0x000001BF2AB... -0.775036 1.089651
4 <rdkit.Chem.rdchem.Mol object at 0x000001BF0E0... 0.346542 -0.820930
y3 y4 y5 y6 y7 y8 y9 ... \
0 -0.427511 0.711080 0.924869 1.096639 1.652419 1.864530 1.112639 ...
1 2.287883 0.424803 -0.721806 -0.398991 -0.414869 0.044067 -0.463779 ...
2 0.167480 0.130518 0.046509 -0.251763 0.650823 0.217649 -0.113869 ...
3 -0.104059 1.107464 1.156964 0.315530 1.703451 1.364106 0.237038 ...
4 -0.902704 0.618991 1.072929 0.839014 -0.128690 0.650888 -1.836564 ...
w5 w6 w7 w8 w9 w10 w11 w12 \
0 1 1 1 1 1 1 1 1
1 1 1 1 1 1 1 1 1
2 1 1 1 1 1 1 1 1
3 1 1 1 1 1 1 1 1
4 1 1 1 1 1 1 1 1
ids split
0 [H]C([H])([H])C([H])([H])[N@@H+]1C([H])([H])[C... train
1 [H]O[C@@]12C([H])([H])O[C@]1([H])[C@]2([H])[C@... train
2 [H]C1([H])C([H])([H])[C@]2([H])[C@@]1([H])[C@]... train
3 [H]OC([H])([H])[C@]1(C([H])([H])C#N)OC1([H])[H] test
4 [H]N1[C@]2([H])C(=O)[C@@]1([H])C([H])([H])OC2(... train
[5 rows x 27 columns]
zmiana nazw kolumn¶
dodaaćlink do data setu + usuniecie kolumn z wagami
# ============================================================
# KOMÓRKA 2: ZMIANA NAZW KOLUMN NA OPISOWE
# ============================================================
# Mapowanie oryginalnych nazw na opisowe
column_mapping = {
# Identyfikatory i features
'ids': 'SMILES',
'X': 'SMILES_raw',
# Właściwości molekularne (y1-y12)
'y1': 'dipole_moment', # μ (Debye) - moment dipolowy
'y2': 'polarizability', # α (Bohr³) - polaryzowalność
'y3': 'HOMO', # εHOMO (eV) - energia HOMO
'y4': 'LUMO', # εLUMO (eV) - energia LUMO
'y5': 'gap', # Δε (eV) - HOMO-LUMO gap
'y6': 'electronic_spatial_extent', # ⟨R²⟩ (Bohr²) - elektroniczny zasięg przestrzenny
'y7': 'zpve', # ZPVE (eV) - zero point vibrational energy
'y8': 'U0', # U₀ (eV) - energia wewnętrzna w 0K
'y9': 'U', # U (eV) - energia wewnętrzna w 298.15K
'y10': 'H', # H (eV) - entalpia w 298.15K
'y11': 'G', # G (eV) - energia swobodna Gibbsa w 298.15K
'y12': 'Cv', # Cv (cal/mol·K) - pojemność cieplna w 298.15K
# Wagi (w1-w12)
'w1': 'weight_dipole',
'w2': 'weight_polarizability',
'w3': 'weight_HOMO',
'w4': 'weight_LUMO',
'w5': 'weight_gap',
'w6': 'weight_R2',
'w7': 'weight_zpve',
'w8': 'weight_U0',
'w9': 'weight_U',
'w10': 'weight_H',
'w11': 'weight_G',
'w12': 'weight_Cv',
# Split (jeśli istnieje)
'split': 'dataset_split'
}
# Zmiana nazw kolumn
df = df.rename(columns=column_mapping)
# ============================================================
# USUNIĘCIE KOLUMN Z WAGAMI (wszystkie = 1)
# ============================================================
# Lista kolumn z wagami do usunięcia
weight_columns = [col for col in df.columns if col.startswith('weight_')]
print("="*70)
print("USUWANIE KOLUMN Z WAGAMI")
print("="*70)
print(f"\n🗑️ Usuwam {len(weight_columns)} kolumn z wagami:")
for col in weight_columns:
print(f" • {col}")
# Usunięcie kolumn
df = df.drop(columns=weight_columns)
print(f"\n✅ Kolumny usunięte! Pozostało {len(df.columns)} kolumn")
print("\n" + "="*70)
print("ZMIENIONO NAZWY KOLUMN")
print("="*70)
print("\n📋 Nowe nazwy kolumn:")
for i, col in enumerate(df.columns, 1):
print(f" {i:2d}. {col}")
print("\n🔍 Pierwsze 3 wiersze po zmianie:")
print(df.head(3))
print("\n📊 Kształt DataFrame:")
print(f" Wiersze: {len(df):,}")
print(f" Kolumny: {len(df.columns)}")
print("\n✅ Nazwy kolumn zostały zmienione!")
print("="*70)
# Opcjonalnie: Wyświetl info o właściwościach
print("\n📖 OPIS WŁAŚCIWOŚCI KWANTOWO-CHEMICZNYCH:")
print("-" * 70)
properties_description = {
'dipole_moment': 'Moment dipolowy [Debye] - miara polarności cząsteczki',
'polarizability': 'Polaryzowalność [Bohr³] - odpowiedź na pole elektryczne',
'HOMO': 'Energia HOMO [eV] - najwyższy zajęty orbital molekularny',
'LUMO': 'Energia LUMO [eV] - najniższy niezajęty orbital molekularny',
'gap': 'HOMO-LUMO gap [eV] - różnica energii, kluczowa dla reaktywności',
'electronic_spatial_extent': 'Elektroniczny zasięg przestrzenny [Bohr²]',
'zpve': 'Energia drgań punktu zerowego [eV]',
'U0': 'Energia wewnętrzna w 0K [eV]',
'U': 'Energia wewnętrzna w 298.15K [eV]',
'H': 'Entalpia w 298.15K [eV]',
'G': 'Energia swobodna Gibbsa w 298.15K [eV]',
'Cv': 'Pojemność cieplna w 298.15K [cal/mol·K]'
}
for prop, desc in properties_description.items():
print(f" • {prop:30s} - {desc}")
print("="*70)
======================================================================
USUWANIE KOLUMN Z WAGAMI
======================================================================
🗑️ Usuwam 12 kolumn z wagami:
• weight_dipole
• weight_polarizability
• weight_HOMO
• weight_LUMO
• weight_gap
• weight_R2
• weight_zpve
• weight_U0
• weight_U
• weight_H
• weight_G
• weight_Cv
✅ Kolumny usunięte! Pozostało 15 kolumn
======================================================================
ZMIENIONO NAZWY KOLUMN
======================================================================
📋 Nowe nazwy kolumn:
1. SMILES_raw
2. dipole_moment
3. polarizability
4. HOMO
5. LUMO
6. gap
7. electronic_spatial_extent
8. zpve
9. U0
10. U
11. H
12. G
13. Cv
14. SMILES
15. dataset_split
🔍 Pierwsze 3 wiersze po zmianie:
SMILES_raw dipole_moment \
0 <rdkit.Chem.rdchem.Mol object at 0x000001BF0ED... -1.291747
1 <rdkit.Chem.rdchem.Mol object at 0x000001BEFC9... -0.644026
2 <rdkit.Chem.rdchem.Mol object at 0x000001BF0E0... -0.459969
polarizability HOMO LUMO gap electronic_spatial_extent \
0 1.802935 -0.427511 0.711080 0.924869 1.096639
1 -0.218808 2.287883 0.424803 -0.721806 -0.398991
2 0.297628 0.167480 0.130518 0.046509 -0.251763
zpve U0 U H G Cv \
0 1.652419 1.864530 1.112639 1.112698 1.112698 1.112566
1 -0.414869 0.044067 -0.463779 -0.463790 -0.463790 -0.463751
2 0.650823 0.217649 -0.113869 -0.113867 -0.113867 -0.113876
SMILES dataset_split
0 [H]C([H])([H])C([H])([H])[N@@H+]1C([H])([H])[C... train
1 [H]O[C@@]12C([H])([H])O[C@]1([H])[C@]2([H])[C@... train
2 [H]C1([H])C([H])([H])[C@]2([H])[C@@]1([H])[C@]... train
📊 Kształt DataFrame:
Wiersze: 1,000
Kolumny: 15
✅ Nazwy kolumn zostały zmienione!
======================================================================
📖 OPIS WŁAŚCIWOŚCI KWANTOWO-CHEMICZNYCH:
----------------------------------------------------------------------
• dipole_moment - Moment dipolowy [Debye] - miara polarności cząsteczki
• polarizability - Polaryzowalność [Bohr³] - odpowiedź na pole elektryczne
• HOMO - Energia HOMO [eV] - najwyższy zajęty orbital molekularny
• LUMO - Energia LUMO [eV] - najniższy niezajęty orbital molekularny
• gap - HOMO-LUMO gap [eV] - różnica energii, kluczowa dla reaktywności
• electronic_spatial_extent - Elektroniczny zasięg przestrzenny [Bohr²]
• zpve - Energia drgań punktu zerowego [eV]
• U0 - Energia wewnętrzna w 0K [eV]
• U - Energia wewnętrzna w 298.15K [eV]
• H - Entalpia w 298.15K [eV]
• G - Energia swobodna Gibbsa w 298.15K [eV]
• Cv - Pojemność cieplna w 298.15K [cal/mol·K]
======================================================================
df.head(10)
| SMILES_raw | dipole_moment | polarizability | HOMO | LUMO | gap | electronic_spatial_extent | zpve | U0 | U | H | G | Cv | SMILES | dataset_split | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | <rdkit.Chem.rdchem.Mol object at 0x000001BF0ED... | -1.291747 | 1.802935 | -0.427511 | 0.711080 | 0.924869 | 1.096639 | 1.652419 | 1.864530 | 1.112639 | 1.112698 | 1.112698 | 1.112566 | [H]C([H])([H])C([H])([H])[N@@H+]1C([H])([H])[C... | train |
| 1 | <rdkit.Chem.rdchem.Mol object at 0x000001BEFC9... | -0.644026 | -0.218808 | 2.287883 | 0.424803 | -0.721806 | -0.398991 | -0.414869 | 0.044067 | -0.463779 | -0.463790 | -0.463790 | -0.463751 | [H]O[C@@]12C([H])([H])O[C@]1([H])[C@]2([H])[C@... | train |
| 2 | <rdkit.Chem.rdchem.Mol object at 0x000001BF0E0... | -0.459969 | 0.297628 | 0.167480 | 0.130518 | 0.046509 | -0.251763 | 0.650823 | 0.217649 | -0.113869 | -0.113867 | -0.113867 | -0.113876 | [H]C1([H])C([H])([H])[C@]2([H])[C@@]1([H])[C@]... | train |
| 3 | <rdkit.Chem.rdchem.Mol object at 0x000001BF2AB... | -0.775036 | 1.089651 | -0.104059 | 1.107464 | 1.156964 | 0.315530 | 1.703451 | 1.364106 | 0.237038 | 0.237068 | 0.237068 | 0.237016 | [H]OC([H])([H])[C@]1(C([H])([H])C#N)OC1([H])[H] | test |
| 4 | <rdkit.Chem.rdchem.Mol object at 0x000001BF0E0... | 0.346542 | -0.820930 | -0.902704 | 0.618991 | 1.072929 | 0.839014 | -0.128690 | 0.650888 | -1.836564 | -1.836554 | -1.836554 | -1.836579 | [H]N1[C@]2([H])C(=O)[C@@]1([H])C([H])([H])OC2(... | train |
| 5 | <rdkit.Chem.rdchem.Mol object at 0x000001BF2A8... | 0.409150 | -2.145600 | -0.758948 | -0.494088 | -0.113557 | -1.271677 | -1.391723 | -1.357742 | -0.503130 | -0.503165 | -0.503165 | -0.503080 | [H][N-]C1OC([H])([H])[C@]12[N@@H+]1C([H])([H])... | valid |
| 6 | <rdkit.Chem.rdchem.Mol object at 0x000001BF2AB... | -0.469601 | 0.294155 | 0.754484 | 0.791158 | 0.412659 | -0.396482 | 0.637064 | 0.321416 | -0.112035 | -0.112037 | -0.112037 | -0.112018 | [H]N([H])C1[NH2+][C@]2([H])C([H])([H])[C@]2(C(... | test |
| 7 | <rdkit.Chem.rdchem.Mol object at 0x000001BEFC5... | 0.327348 | -0.128489 | -0.611199 | -1.142716 | -0.835853 | -0.335116 | -1.578459 | -1.002448 | -0.027435 | -0.027466 | -0.027466 | -0.027395 | [H]C1([H])[N@@H+]2[C@@]3([H])[C@]2([H])[C@@]2(... | train |
| 8 | <rdkit.Chem.rdchem.Mol object at 0x000001BF2A2... | 1.261589 | -2.541611 | -1.797187 | -0.974553 | -0.073540 | -1.369651 | -1.949059 | -2.116866 | 0.891255 | 0.891221 | 0.891221 | 0.891306 | [H]C([H])([H])O[C@@]1(C([H])([H])[H])C([H])([H... | train |
| 9 | <rdkit.Chem.rdchem.Mol object at 0x000001BF0E1... | 1.653432 | -0.118068 | 0.462979 | -1.759313 | -1.990326 | -0.729647 | -1.155817 | -1.150208 | -0.435730 | -0.435769 | -0.435769 | -0.435686 | [H]C(=O)c1c([H])c(N([H])[H])nn1[H] | train |
zmana 1 kolumny na licaby¶
# ============================================================
# KOMÓRKA 3: USUNIĘCIE ZBĘDNEJ KOLUMNY SMILES_raw
# ============================================================
print("="*70)
print("USUWANIE ZBĘDNEJ KOLUMNY")
print("="*70)
# Sprawdzenie czy kolumna istnieje
if 'SMILES_raw' in df.columns:
print("\n🗑️ Usuwam kolumnę: SMILES_raw (duplikat SMILES)")
df = df.drop(columns=['SMILES_raw'])
print("✅ Kolumna usunięta!")
else:
print("\n⚠️ Kolumna SMILES_raw nie istnieje")
print(f"\n📊 Aktualne kolumny ({len(df.columns)}):")
for i, col in enumerate(df.columns, 1):
print(f" {i:2d}. {col}")
print(f"\n📈 Kształt DataFrame: {df.shape}")
print("="*70)
====================================================================== USUWANIE ZBĘDNEJ KOLUMNY ====================================================================== 🗑️ Usuwam kolumnę: SMILES_raw (duplikat SMILES) ✅ Kolumna usunięta! 📊 Aktualne kolumny (14): 1. dipole_moment 2. polarizability 3. HOMO 4. LUMO 5. gap 6. electronic_spatial_extent 7. zpve 8. U0 9. U 10. H 11. G 12. Cv 13. SMILES 14. dataset_split 📈 Kształt DataFrame: (1000, 14) ======================================================================
df.head(10)
| dipole_moment | polarizability | HOMO | LUMO | gap | electronic_spatial_extent | zpve | U0 | U | H | G | Cv | SMILES | dataset_split | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | -1.291747 | 1.802935 | -0.427511 | 0.711080 | 0.924869 | 1.096639 | 1.652419 | 1.864530 | 1.112639 | 1.112698 | 1.112698 | 1.112566 | [H]C([H])([H])C([H])([H])[N@@H+]1C([H])([H])[C... | train |
| 1 | -0.644026 | -0.218808 | 2.287883 | 0.424803 | -0.721806 | -0.398991 | -0.414869 | 0.044067 | -0.463779 | -0.463790 | -0.463790 | -0.463751 | [H]O[C@@]12C([H])([H])O[C@]1([H])[C@]2([H])[C@... | train |
| 2 | -0.459969 | 0.297628 | 0.167480 | 0.130518 | 0.046509 | -0.251763 | 0.650823 | 0.217649 | -0.113869 | -0.113867 | -0.113867 | -0.113876 | [H]C1([H])C([H])([H])[C@]2([H])[C@@]1([H])[C@]... | train |
| 3 | -0.775036 | 1.089651 | -0.104059 | 1.107464 | 1.156964 | 0.315530 | 1.703451 | 1.364106 | 0.237038 | 0.237068 | 0.237068 | 0.237016 | [H]OC([H])([H])[C@]1(C([H])([H])C#N)OC1([H])[H] | test |
| 4 | 0.346542 | -0.820930 | -0.902704 | 0.618991 | 1.072929 | 0.839014 | -0.128690 | 0.650888 | -1.836564 | -1.836554 | -1.836554 | -1.836579 | [H]N1[C@]2([H])C(=O)[C@@]1([H])C([H])([H])OC2(... | train |
| 5 | 0.409150 | -2.145600 | -0.758948 | -0.494088 | -0.113557 | -1.271677 | -1.391723 | -1.357742 | -0.503130 | -0.503165 | -0.503165 | -0.503080 | [H][N-]C1OC([H])([H])[C@]12[N@@H+]1C([H])([H])... | valid |
| 6 | -0.469601 | 0.294155 | 0.754484 | 0.791158 | 0.412659 | -0.396482 | 0.637064 | 0.321416 | -0.112035 | -0.112037 | -0.112037 | -0.112018 | [H]N([H])C1[NH2+][C@]2([H])C([H])([H])[C@]2(C(... | test |
| 7 | 0.327348 | -0.128489 | -0.611199 | -1.142716 | -0.835853 | -0.335116 | -1.578459 | -1.002448 | -0.027435 | -0.027466 | -0.027466 | -0.027395 | [H]C1([H])[N@@H+]2[C@@]3([H])[C@]2([H])[C@@]2(... | train |
| 8 | 1.261589 | -2.541611 | -1.797187 | -0.974553 | -0.073540 | -1.369651 | -1.949059 | -2.116866 | 0.891255 | 0.891221 | 0.891221 | 0.891306 | [H]C([H])([H])O[C@@]1(C([H])([H])[H])C([H])([H... | train |
| 9 | 1.653432 | -0.118068 | 0.462979 | -1.759313 | -1.990326 | -0.729647 | -1.155817 | -1.150208 | -0.435730 | -0.435769 | -0.435769 | -0.435686 | [H]C(=O)c1c([H])c(N([H])[H])nn1[H] | train |
Wizualizacja¶
# ============================================================
# KOMÓRKA 4: WIZUALIZACJA STRUKTUR MOLEKULARNYCH
# ============================================================
from rdkit import Chem
from rdkit.Chem import Draw
import matplotlib.pyplot as plt
from matplotlib.gridspec import GridSpec
print("="*70)
print("WIZUALIZACJA PRZYKŁADOWYCH CZĄSTECZEK")
print("="*70)
# ============================================================
# Wybór przykładowych cząsteczek
# ============================================================
# Losowe próbki z datasetu
n_samples = 12
sample_molecules = df.sample(n=n_samples, random_state=42)
print(f"\n🔬 Wybieram {n_samples} losowych cząsteczek do wizualizacji...\n")
# ============================================================
# Przygotowanie molekuł RDKit
# ============================================================
molecules = []
labels = []
for idx, row in sample_molecules.iterrows():
smiles = row['SMILES']
mol = Chem.MolFromSmiles(smiles)
if mol is not None:
# Informacje o cząsteczce
gap = row['gap'] if 'gap' in row else None
homo = row['HOMO'] if 'HOMO' in row else None
# Label z kluczowymi informacjami
label = f"ID: {idx}\n"
label += f"SMILES: {smiles[:20]}{'...' if len(smiles) > 20 else ''}\n"
if gap is not None:
label += f"Gap: {gap:.3f} eV"
molecules.append(mol)
labels.append(label)
print(f" ✓ Cząsteczka {len(molecules):2d}: {smiles[:40]}{'...' if len(smiles) > 40 else ''}")
print(f"\n✅ Przygotowano {len(molecules)} struktur do wizualizacji")
# ============================================================
# Wizualizacja - Grid 3x4
# ============================================================
print("\n🎨 Generowanie wizualizacji...\n")
fig = plt.figure(figsize=(16, 12))
gs = GridSpec(3, 4, figure=fig, hspace=0.4, wspace=0.3)
for i, (mol, label) in enumerate(zip(molecules, labels)):
row = i // 4
col = i % 4
ax = fig.add_subplot(gs[row, col])
# Rysowanie cząsteczki
img = Draw.MolToImage(mol, size=(400, 400))
ax.imshow(img)
ax.axis('off')
ax.set_title(label, fontsize=9, pad=10)
plt.suptitle('Przykładowe cząsteczki z datasetu QM9',
fontsize=16, fontweight='bold', y=0.98)
plt.tight_layout()
plt.savefig('qm9_example_molecules.png', dpi=300, bbox_inches='tight')
print("✅ Zapisano: qm9_example_molecules.png")
plt.show()
# ============================================================
# Wizualizacja cząsteczek o ekstremalnych właściwościach
# ============================================================
print("\n" + "="*70)
print("CZĄSTECZKI O EKSTREMALNYCH WŁAŚCIWOŚCIACH")
print("="*70)
# Najwyższy i najniższy gap
idx_max_gap = df['gap'].idxmax()
idx_min_gap = df['gap'].idxmin()
# Najwyższy i najniższy moment dipolowy
idx_max_dipole = df['dipole_moment'].idxmax()
idx_min_dipole = df['dipole_moment'].idxmin()
extreme_indices = [idx_max_gap, idx_min_gap, idx_max_dipole, idx_min_dipole]
extreme_labels_text = [
f"Najwyższy gap: {df.loc[idx_max_gap, 'gap']:.3f} eV",
f"Najniższy gap: {df.loc[idx_min_gap, 'gap']:.3f} eV",
f"Najwyższy dipole: {df.loc[idx_max_dipole, 'dipole_moment']:.3f} D",
f"Najniższy dipole: {df.loc[idx_min_dipole, 'dipole_moment']:.3f} D"
]
print("\n🔍 Ekstremalne cząsteczki:")
for i, (idx, text) in enumerate(zip(extreme_indices, extreme_labels_text), 1):
smiles = df.loc[idx, 'SMILES']
print(f" {i}. {text}")
print(f" SMILES: {smiles}")
# Wizualizacja ekstremalnych cząsteczek
fig, axes = plt.subplots(2, 2, figsize=(12, 12))
axes = axes.flatten()
for i, (idx, label_text) in enumerate(zip(extreme_indices, extreme_labels_text)):
smiles = df.loc[idx, 'SMILES']
mol = Chem.MolFromSmiles(smiles)
if mol is not None:
img = Draw.MolToImage(mol, size=(500, 500))
axes[i].imshow(img)
axes[i].axis('off')
axes[i].set_title(f"{label_text}\nSMILES: {smiles[:30]}{'...' if len(smiles) > 30 else ''}",
fontsize=10, fontweight='bold', pad=10)
plt.suptitle('Cząsteczki o ekstremalnych właściwościach',
fontsize=14, fontweight='bold', y=0.98)
plt.tight_layout()
plt.savefig('qm9_extreme_molecules.png', dpi=300, bbox_inches='tight')
print("\n✅ Zapisano: qm9_extreme_molecules.png")
plt.show()
# ============================================================
# Przykłady według rozmiaru
# ============================================================
print("\n" + "="*70)
print("CZĄSTECZKI WEDŁUG ROZMIARU")
print("="*70)
# Dodanie kolumny z liczbą atomów
df['num_atoms'] = df['SMILES'].apply(lambda x: Chem.MolFromSmiles(x).GetNumAtoms()
if Chem.MolFromSmiles(x) is not None else 0)
# Najmniejsze i największe
idx_smallest = df['num_atoms'].idxmin()
idx_largest = df['num_atoms'].idxmax()
# Średnie
median_atoms = df['num_atoms'].median()
idx_medium = df.iloc[(df['num_atoms'] - median_atoms).abs().argsort()[:1]].index[0]
size_indices = [idx_smallest, idx_medium, idx_largest]
size_labels = [
f"Najmniejsza ({df.loc[idx_smallest, 'num_atoms']} atomów)",
f"Średnia ({df.loc[idx_medium, 'num_atoms']} atomów)",
f"Największa ({df.loc[idx_largest, 'num_atoms']} atomów)"
]
print("\n📏 Cząsteczki według rozmiaru:")
for label, idx in zip(size_labels, size_indices):
smiles = df.loc[idx, 'SMILES']
print(f" • {label}: {smiles}")
# Wizualizacja
fig, axes = plt.subplots(1, 3, figsize=(15, 5))
for i, (idx, label) in enumerate(zip(size_indices, size_labels)):
smiles = df.loc[idx, 'SMILES']
mol = Chem.MolFromSmiles(smiles)
if mol is not None:
img = Draw.MolToImage(mol, size=(500, 500))
axes[i].imshow(img)
axes[i].axis('off')
axes[i].set_title(f"{label}\n{smiles}", fontsize=10, fontweight='bold')
plt.suptitle('Porównanie rozmiaru cząsteczek', fontsize=14, fontweight='bold')
plt.tight_layout()
plt.savefig('qm9_size_comparison.png', dpi=300, bbox_inches='tight')
print("\n✅ Zapisano: qm9_size_comparison.png")
plt.show()
print("\n" + "="*70)
print("✅ WIZUALIZACJE ZAKOŃCZONE!")
print(" Wygenerowano 3 pliki PNG z strukturami molekularnymi")
print("="*70)
====================================================================== WIZUALIZACJA PRZYKŁADOWYCH CZĄSTECZEK ====================================================================== 🔬 Wybieram 12 losowych cząsteczek do wizualizacji... ✓ Cząsteczka 1: [H]Oc1c([H])c(N([H])[H])n([H])c1C([H])([... ✓ Cząsteczka 2: [H]O[C@]1([H])C([H])([H])[N@@H+]2C([H])(... ✓ Cząsteczka 3: [H]C(=O)[C@]1(C([H])([H])[H])C([H])([H])... ✓ Cząsteczka 4: [H]OC([H])([H])[C@](O[H])(C([H])([H])[H]... ✓ Cząsteczka 5: [H]OC([H])([H])[C@]12[C@@]3([H])[C@]4([H... ✓ Cząsteczka 6: [H][N-]C1O[C@@]2([H])[C@@]([H])(O1)[C@@]... ✓ Cząsteczka 7: [H]C([H])([H])OC(=O)[C@@]([H])(OC([H])([... ✓ Cząsteczka 8: [H]C1([H])O[C@@]12[C@@]1([H])[N@H+]3[C@]... ✓ Cząsteczka 9: [H]OC([H])([H])C([H])([H])N(C([H])=O)C([... ✓ Cząsteczka 10: [H]C1=C([H])[C@@]([H])(C([H])([H])[H])OC... ✓ Cząsteczka 11: [H]C#C[C@@]1([C@]([H])(O[H])C([H])([H])O... ✓ Cząsteczka 12: [H]Oc1nc(N([H])C([H])([H])[H])oc1C([H])(... ✅ Przygotowano 12 struktur do wizualizacji 🎨 Generowanie wizualizacji...
C:\Users\slast\AppData\Local\Temp\ipykernel_6436\3317227650.py:77: UserWarning: This figure includes Axes that are not compatible with tight_layout, so results might be incorrect. plt.tight_layout()
✅ Zapisano: qm9_example_molecules.png
======================================================================
CZĄSTECZKI O EKSTREMALNYCH WŁAŚCIWOŚCIACH
======================================================================
🔍 Ekstremalne cząsteczki:
1. Najwyższy gap: 4.752 eV
SMILES: [H]C#C[C@@]1(C([H])([H])C([H])=O)N([H])[C@]1([H])C([H])([H])[H]
2. Najniższy gap: -2.537 eV
SMILES: [H]C1=C([H])c2c(oc([H])c2[H])C1=O
3. Najwyższy dipole: 4.349 D
SMILES: [H]C1=C([H])c2c(oc([H])c2[H])C1=O
4. Najniższy dipole: -1.694 D
SMILES: [H]C#C[C@@]1(C([H])([H])C([H])=O)N([H])[C@]1([H])C([H])([H])[H]
✅ Zapisano: qm9_extreme_molecules.png
====================================================================== CZĄSTECZKI WEDŁUG ROZMIARU ====================================================================== 📏 Cząsteczki według rozmiaru: • Najmniejsza (4 atomów): [H]N([H])C(=O)C([H])([H])[H] • Średnia (9 atomów): [H]O[C@@]12C([H])([H])O[C@]1([H])[C@]2([H])[C@]([H])(O[H])C([H])([H])[H] • Największa (11 atomów): [H]/N=C(/C(=O)N(/C([H])=N/[H])C([H])([H])[H])N([H])[H] ✅ Zapisano: qm9_size_comparison.png
====================================================================== ✅ WIZUALIZACJE ZAKOŃCZONE! Wygenerowano 3 pliki PNG z strukturami molekularnymi ======================================================================
brakujace wartosaci i duplikaty¶
# ============================================================
# SPRAWDZANIE BRAKUJĄCYCH WARTOŚCI
# ============================================================
missing_counts = df.isnull().sum()
total_missing = missing_counts.sum()
if total_missing == 0:
print("✔️ Brak brakujących wartości w zbiorze danych.")
else:
print(f"⚠️ Znaleziono {total_missing} brakujących wartości w zbiorze danych:")
print(missing_counts[missing_counts > 0])
# ============================================================
# SPRAWDZANIE DUPLIKATÓW
# ============================================================
duplicate_count = df.duplicated().sum()
if duplicate_count == 0:
print("✔️ Brak duplikatów w zbiorze danych.")
else:
print(f"⚠️ Znaleziono {duplicate_count} duplikatów w zbiorze danych.")
✔️ Brak brakujących wartości w zbiorze danych. ✔️ Brak duplikatów w zbiorze danych.
teraz czas na sprawdzenie czy mamy jkaieś odstające wartosci i wyrysoanie box plotów¶
import matplotlib.pyplot as plt
import seaborn as sns
# ============================================================
# WYBIERAMY TYLKO KOLUMNY NUMERYCZNE
# ============================================================
numeric_cols = df.select_dtypes(include=['float64', 'int64']).columns
print("📊 Sprawdzenie wartości odstających w kolumnach numerycznych:\n")
for col in numeric_cols:
Q1 = df[col].quantile(0.25)
Q3 = df[col].quantile(0.75)
IQR = Q3 - Q1
lower_bound = Q1 - 1.5 * IQR
upper_bound = Q3 + 1.5 * IQR
outliers = df[(df[col] < lower_bound) | (df[col] > upper_bound)]
print(f"{col}: {len(outliers)} wartości odstających")
# ============================================================
# BOX PLOTY DLA KOLUMN NUMERYCZNYCH
# ============================================================
plt.figure(figsize=(15, len(numeric_cols)*1.5))
for i, col in enumerate(numeric_cols, 1):
plt.subplot(len(numeric_cols), 1, i)
sns.boxplot(x=df[col], color='skyblue')
plt.title(f"Boxplot: {col}")
plt.tight_layout()
plt.show()
📊 Sprawdzenie wartości odstających w kolumnach numerycznych: dipole_moment: 10 wartości odstających polarizability: 33 wartości odstających HOMO: 68 wartości odstających LUMO: 0 wartości odstających gap: 2 wartości odstających electronic_spatial_extent: 70 wartości odstających zpve: 3 wartości odstających U0: 17 wartości odstających U: 25 wartości odstających H: 25 wartości odstających G: 25 wartości odstających Cv: 25 wartości odstających num_atoms: 195 wartości odstających
wzory kilku czasteczek o wartosciach odstajacych¶
print(df.columns)
Index(['dipole_moment', 'polarizability', 'HOMO', 'LUMO', 'gap',
'electronic_spatial_extent', 'zpve', 'U0', 'U', 'H', 'G', 'Cv',
'SMILES', 'dataset_split', 'num_atoms'],
dtype='object')
# ============================================================
# WYŚWIETLANIE CZĄSTECZEK Z WARTOŚCIAMI ODSTAJĄCYMI
# ============================================================
from rdkit import Chem
from rdkit.Chem import Draw
# Nazwa kolumny SMILES w Twoim df
smiles_col = 'SMILES'
# Sprawdzenie, które kolumny numeryczne mają wartości odstające
numeric_cols = df.select_dtypes(include=['float64', 'int64']).columns
# Zbieramy wszystkie indeksy z wartościami odstającymi
outlier_indices = set()
for col in numeric_cols:
Q1 = df[col].quantile(0.25)
Q3 = df[col].quantile(0.75)
IQR = Q3 - Q1
lower_bound = Q1 - 1.5 * IQR
upper_bound = Q3 + 1.5 * IQR
outliers = df[(df[col] < lower_bound) | (df[col] > upper_bound)]
outlier_indices.update(outliers.index.tolist())
# Wybieramy np. pierwsze 5 cząsteczek z wartościami odstającymi
sample_outliers = df.loc[list(outlier_indices)[:5], smiles_col]
# Konwersja SMILES na obiekty RDKit
mols = [Chem.MolFromSmiles(s) for s in sample_outliers]
# Wyświetlenie cząsteczek
Draw.MolsToGridImage(
mols,
molsPerRow=5,
subImgSize=(200,200),
legends=[f"{i+1}" for i in range(len(mols))]
)
# ============================================================
# KOMÓRKA 5: ANALIZA KORELACJI MIĘDZY WŁAŚCIWOŚCIAMI
# ============================================================
import numpy as np
import seaborn as sns
import matplotlib.pyplot as plt
print("="*70)
print("ANALIZA KORELACJI - WŁAŚCIWOŚCI QM9")
print("="*70)
# ============================================================
# Przygotowanie danych numerycznych
# ============================================================
# Wybór tylko kolumn numerycznych (bez SMILES i split)
numeric_cols = df.select_dtypes(include=[np.number]).columns.tolist()
# Usunięcie kolumny num_atoms jeśli istnieje (dodana w poprzedniej komórce)
if 'num_atoms' in numeric_cols:
numeric_cols.remove('num_atoms')
print(f"\n📊 Analizowane właściwości ({len(numeric_cols)}):")
for i, col in enumerate(numeric_cols, 1):
print(f" {i:2d}. {col}")
df_numeric = df[numeric_cols].copy()
# ============================================================
# 1. MACIERZ KORELACJI - PEŁNA
# ============================================================
print("\n" + "="*70)
print("1. PEŁNA MACIERZ KORELACJI PEARSONA")
print("="*70)
# Obliczenie korelacji
correlation_matrix = df_numeric.corr(method='pearson')
print(f"\n✓ Macierz korelacji obliczona: {correlation_matrix.shape}")
# Wizualizacja - pełna macierz
plt.figure(figsize=(14, 12))
mask = np.triu(np.ones_like(correlation_matrix, dtype=bool), k=1)
sns.heatmap(correlation_matrix,
mask=mask,
annot=True,
fmt='.2f',
cmap='RdBu_r',
center=0,
square=True,
linewidths=0.5,
cbar_kws={"shrink": 0.8, "label": "Korelacja Pearsona"},
vmin=-1,
vmax=1)
plt.title('Macierz korelacji właściwości kwantowo-chemicznych\n(Korelacja Pearsona)',
fontsize=14, fontweight='bold', pad=20)
plt.xticks(rotation=45, ha='right')
plt.yticks(rotation=0)
plt.tight_layout()
plt.savefig('qm9_correlation_full.png', dpi=300, bbox_inches='tight')
print("✅ Zapisano: qm9_correlation_full.png")
plt.show()
# ============================================================
# 2. NAJSILNIEJSZE KORELACJE
# ============================================================
print("\n" + "="*70)
print("2. TOP 15 NAJSILNIEJSZYCH KORELACJI")
print("="*70)
# Ekstrakcja par korelacji (bez diagonali i duplikatów)
correlations_list = []
for i in range(len(correlation_matrix.columns)):
for j in range(i+1, len(correlation_matrix.columns)):
correlations_list.append({
'Feature 1': correlation_matrix.columns[i],
'Feature 2': correlation_matrix.columns[j],
'Correlation': correlation_matrix.iloc[i, j]
})
df_correlations = pd.DataFrame(correlations_list)
df_correlations['Abs_Correlation'] = df_correlations['Correlation'].abs()
df_correlations = df_correlations.sort_values('Abs_Correlation', ascending=False)
print("\n🔝 TOP 15 najsilniejszych korelacji (wartość bezwzględna):")
print("-" * 70)
for idx, row in df_correlations.head(15).iterrows():
print(f" {row['Feature 1']:25s} ↔ {row['Feature 2']:25s} r = {row['Correlation']:+.4f}")
# Wizualizacja top korelacji
top_n = 15
top_correlations = df_correlations.head(top_n).copy()
top_correlations['Pair'] = (top_correlations['Feature 1'] + '\n↔\n' +
top_correlations['Feature 2'])
plt.figure(figsize=(12, 8))
colors = ['#d73027' if x < 0 else '#4575b4' for x in top_correlations['Correlation']]
bars = plt.barh(range(len(top_correlations)),
top_correlations['Correlation'],
color=colors,
edgecolor='black',
linewidth=0.7)
plt.yticks(range(len(top_correlations)), top_correlations['Pair'], fontsize=9)
plt.xlabel('Współczynnik korelacji Pearsona', fontsize=12, fontweight='bold')
plt.title('Top 15 najsilniejszych korelacji między właściwościami',
fontsize=14, fontweight='bold', pad=20)
plt.axvline(x=0, color='black', linewidth=0.8, linestyle='-')
plt.grid(axis='x', alpha=0.3, linestyle='--')
# Dodanie wartości na słupkach
for i, (bar, val) in enumerate(zip(bars, top_correlations['Correlation'])):
plt.text(val + 0.02 if val > 0 else val - 0.02,
i,
f'{val:.3f}',
va='center',
ha='left' if val > 0 else 'right',
fontweight='bold',
fontsize=9)
plt.tight_layout()
plt.savefig('qm9_top_correlations.png', dpi=300, bbox_inches='tight')
print("\n✅ Zapisano: qm9_top_correlations.png")
plt.show()
# ============================================================
# 3. KORELACJA HOMO-LUMO-GAP (kluczowe właściwości)
# ============================================================
print("\n" + "="*70)
print("3. ANALIZA HOMO-LUMO-GAP")
print("="*70)
# Korelacje z gap
gap_correlations = correlation_matrix['gap'].sort_values(ascending=False)
print("\n📊 Korelacje z HOMO-LUMO gap:")
print("-" * 70)
for prop, corr in gap_correlations.items():
if prop != 'gap':
print(f" {prop:30s} r = {corr:+.4f}")
# Scatter plots - HOMO, LUMO, gap
fig, axes = plt.subplots(1, 3, figsize=(18, 5))
# Sample dla szybszej wizualizacji
df_sample = df_numeric.sample(n=min(5000, len(df_numeric)), random_state=42)
# HOMO vs LUMO
axes[0].scatter(df_sample['HOMO'], df_sample['LUMO'],
alpha=0.3, s=10, c='#1f77b4', edgecolors='none')
axes[0].set_xlabel('HOMO (eV)', fontsize=11, fontweight='bold')
axes[0].set_ylabel('LUMO (eV)', fontsize=11, fontweight='bold')
axes[0].set_title(f'HOMO vs LUMO\nr = {correlation_matrix.loc["HOMO", "LUMO"]:.3f}',
fontsize=12, fontweight='bold')
axes[0].grid(alpha=0.3)
# HOMO vs gap
axes[1].scatter(df_sample['HOMO'], df_sample['gap'],
alpha=0.3, s=10, c='#ff7f0e', edgecolors='none')
axes[1].set_xlabel('HOMO (eV)', fontsize=11, fontweight='bold')
axes[1].set_ylabel('gap (eV)', fontsize=11, fontweight='bold')
axes[1].set_title(f'HOMO vs gap\nr = {correlation_matrix.loc["HOMO", "gap"]:.3f}',
fontsize=12, fontweight='bold')
axes[1].grid(alpha=0.3)
# LUMO vs gap
axes[2].scatter(df_sample['LUMO'], df_sample['gap'],
alpha=0.3, s=10, c='#2ca02c', edgecolors='none')
axes[2].set_xlabel('LUMO (eV)', fontsize=11, fontweight='bold')
axes[2].set_ylabel('gap (eV)', fontsize=11, fontweight='bold')
axes[2].set_title(f'LUMO vs gap\nr = {correlation_matrix.loc["LUMO", "gap"]:.3f}',
fontsize=12, fontweight='bold')
axes[2].grid(alpha=0.3)
plt.suptitle('Zależności między kluczowymi właściwościami elektronowymi',
fontsize=14, fontweight='bold', y=1.02)
plt.tight_layout()
plt.savefig('qm9_homo_lumo_gap.png', dpi=300, bbox_inches='tight')
print("\n✅ Zapisano: qm9_homo_lumo_gap.png")
plt.show()
# ============================================================
# 4. PAIRPLOT - WYBRANE WŁAŚCIWOŚCI
# ============================================================
print("\n" + "="*70)
print("4. PAIRPLOT - KLUCZOWE WŁAŚCIWOŚCI")
print("="*70)
# Wybór najważniejszych właściwości
key_properties = ['HOMO', 'LUMO', 'gap', 'dipole_moment', 'polarizability']
print(f"\n📈 Tworzenie pairplot dla {len(key_properties)} właściwości...")
print(f" Używam {min(2000, len(df))} próbek dla szybszej wizualizacji")
df_pairplot = df[key_properties].sample(n=min(2000, len(df)), random_state=42)
# Pairplot
pairplot_fig = sns.pairplot(df_pairplot,
diag_kind='kde',
plot_kws={'alpha': 0.4, 's': 15, 'edgecolor': 'none'},
diag_kws={'linewidth': 2})
pairplot_fig.fig.suptitle('Pairplot kluczowych właściwości kwantowo-chemicznych',
fontsize=14, fontweight='bold', y=1.01)
plt.savefig('qm9_pairplot.png', dpi=300, bbox_inches='tight')
print("✅ Zapisano: qm9_pairplot.png")
plt.show()
# ============================================================
# 5. KORELACJA SPEARMANA (dla sprawdzenia nieliniowych zależności)
# ============================================================
print("\n" + "="*70)
print("5. PORÓWNANIE KORELACJI PEARSONA vs SPEARMANA")
print("="*70)
# Korelacja Spearmana
spearman_matrix = df_numeric.corr(method='spearman')
# Różnica między Pearson a Spearman
diff_matrix = np.abs(correlation_matrix - spearman_matrix)
print("\n📊 Największe różnice między Pearson a Spearman (nieliniowości):")
print("-" * 70)
# Znalezienie największych różnic
max_diffs = []
for i in range(len(diff_matrix.columns)):
for j in range(i+1, len(diff_matrix.columns)):
max_diffs.append({
'Feature 1': diff_matrix.columns[i],
'Feature 2': diff_matrix.columns[j],
'Pearson': correlation_matrix.iloc[i, j],
'Spearman': spearman_matrix.iloc[i, j],
'Difference': diff_matrix.iloc[i, j]
})
df_diffs = pd.DataFrame(max_diffs).sort_values('Difference', ascending=False)
for idx, row in df_diffs.head(10).iterrows():
print(f" {row['Feature 1']:20s} ↔ {row['Feature 2']:20s}")
print(f" Pearson: {row['Pearson']:+.4f} | Spearman: {row['Spearman']:+.4f} | Diff: {row['Difference']:.4f}")
# ============================================================
# 6. PODSUMOWANIE STATYSTYCZNE
# ============================================================
print("\n" + "="*70)
print("6. PODSUMOWANIE ANALIZY KORELACJI")
print("="*70)
print(f"""
📊 STATYSTYKI KORELACJI:
• Liczba analizowanych właściwości: {len(numeric_cols)}
• Liczba par korelacji: {len(correlations_list)}
• Najsilniejsza korelacja dodatnia:
{df_correlations.iloc[0]['Feature 1']} ↔ {df_correlations.iloc[0]['Feature 2']}
r = {df_correlations.iloc[0]['Correlation']:.4f}
• Najsilniejsza korelacja ujemna:
{df_correlations[df_correlations['Correlation'] < 0].iloc[0]['Feature 1']} ↔ {df_correlations[df_correlations['Correlation'] < 0].iloc[0]['Feature 2']}
r = {df_correlations[df_correlations['Correlation'] < 0].iloc[0]['Correlation']:.4f}
• Silne korelacje (|r| > 0.7): {len(df_correlations[df_correlations['Abs_Correlation'] > 0.7])}
• Średnie korelacje (0.3 < |r| < 0.7): {len(df_correlations[(df_correlations['Abs_Correlation'] > 0.3) & (df_correlations['Abs_Correlation'] <= 0.7)])}
• Słabe korelacje (|r| < 0.3): {len(df_correlations[df_correlations['Abs_Correlation'] <= 0.3])}
""")
print("="*70)
print("✅ ANALIZA KORELACJI ZAKOŃCZONA!")
print(" Wygenerowano 4 pliki PNG z wizualizacjami")
print("="*70)
====================================================================== ANALIZA KORELACJI - WŁAŚCIWOŚCI QM9 ====================================================================== 📊 Analizowane właściwości (12): 1. dipole_moment 2. polarizability 3. HOMO 4. LUMO 5. gap 6. electronic_spatial_extent 7. zpve 8. U0 9. U 10. H 11. G 12. Cv ====================================================================== 1. PEŁNA MACIERZ KORELACJI PEARSONA ====================================================================== ✓ Macierz korelacji obliczona: (12, 12) ✅ Zapisano: qm9_correlation_full.png
====================================================================== 2. TOP 15 NAJSILNIEJSZYCH KORELACJI ====================================================================== 🔝 TOP 15 najsilniejszych korelacji (wartość bezwzględna): ---------------------------------------------------------------------- H ↔ G r = +1.0000 U ↔ H r = +1.0000 U ↔ G r = +1.0000 U ↔ Cv r = +1.0000 H ↔ Cv r = +1.0000 G ↔ Cv r = +1.0000 LUMO ↔ gap r = +0.8670 electronic_spatial_extent ↔ U0 r = +0.8288 polarizability ↔ U0 r = +0.7982 zpve ↔ U0 r = +0.7758 polarizability ↔ zpve r = +0.7655 polarizability ↔ electronic_spatial_extent r = +0.7170 LUMO ↔ zpve r = +0.6705 gap ↔ zpve r = +0.5705 electronic_spatial_extent ↔ zpve r = +0.5284 ✅ Zapisano: qm9_top_correlations.png
====================================================================== 3. ANALIZA HOMO-LUMO-GAP ====================================================================== 📊 Korelacje z HOMO-LUMO gap: ---------------------------------------------------------------------- LUMO r = +0.8670 zpve r = +0.5705 G r = +0.3873 H r = +0.3873 U r = +0.3873 Cv r = +0.3873 U0 r = +0.2404 polarizability r = +0.1883 electronic_spatial_extent r = +0.0438 HOMO r = -0.2878 dipole_moment r = -0.3407 ✅ Zapisano: qm9_homo_lumo_gap.png
====================================================================== 4. PAIRPLOT - KLUCZOWE WŁAŚCIWOŚCI ====================================================================== 📈 Tworzenie pairplot dla 5 właściwości... Używam 1000 próbek dla szybszej wizualizacji ✅ Zapisano: qm9_pairplot.png
======================================================================
5. PORÓWNANIE KORELACJI PEARSONA vs SPEARMANA
======================================================================
📊 Największe różnice między Pearson a Spearman (nieliniowości):
----------------------------------------------------------------------
polarizability ↔ H
Pearson: -0.1031 | Spearman: +0.1497 | Diff: 0.2528
polarizability ↔ G
Pearson: -0.1031 | Spearman: +0.1497 | Diff: 0.2528
polarizability ↔ U
Pearson: -0.1031 | Spearman: +0.1496 | Diff: 0.2527
polarizability ↔ Cv
Pearson: -0.1031 | Spearman: +0.1495 | Diff: 0.2526
U0 ↔ H
Pearson: -0.3127 | Spearman: -0.1948 | Diff: 0.1179
U0 ↔ G
Pearson: -0.3127 | Spearman: -0.1948 | Diff: 0.1179
U0 ↔ U
Pearson: -0.3127 | Spearman: -0.1950 | Diff: 0.1177
U0 ↔ Cv
Pearson: -0.3128 | Spearman: -0.1953 | Diff: 0.1175
polarizability ↔ gap
Pearson: +0.1883 | Spearman: +0.2940 | Diff: 0.1057
HOMO ↔ U0
Pearson: +0.1213 | Spearman: +0.0184 | Diff: 0.1028
======================================================================
6. PODSUMOWANIE ANALIZY KORELACJI
======================================================================
📊 STATYSTYKI KORELACJI:
• Liczba analizowanych właściwości: 12
• Liczba par korelacji: 66
• Najsilniejsza korelacja dodatnia:
H ↔ G
r = 1.0000
• Najsilniejsza korelacja ujemna:
dipole_moment ↔ LUMO
r = -0.3941
• Silne korelacje (|r| > 0.7): 12
• Średnie korelacje (0.3 < |r| < 0.7): 24
• Słabe korelacje (|r| < 0.3): 30
======================================================================
✅ ANALIZA KORELACJI ZAKOŃCZONA!
Wygenerowano 4 pliki PNG z wizualizacjami
======================================================================
🔝 TOP 15 NAJSILNIEJSZYCH KORELACJI¶
Poniżej przedstawiono 15 par właściwości z najsilniejszymi korelacjami w datasecie QM9:
Korelacje perfekcyjne (r = 1.00)¶
Te właściwości są ze sobą matematycznie związane poprzez definicje termodynamiczne:
- H ↔ G → r = +1.0000 (Entalpia vs Energia Gibbsa)
- U ↔ H → r = +1.0000 (Energia wewnętrzna vs Entalpia)
- U ↔ G → r = +1.0000 (Energia wewnętrzna vs Energia Gibbsa)
- U ↔ Cv → r = +1.0000 (Energia wewnętrzna vs Pojemność cieplna)
- H ↔ Cv → r = +1.0000 (Entalpia vs Pojemność cieplna)
- G ↔ Cv → r = +1.0000 (Energia Gibbsa vs Pojemność cieplna)
⚠️ Uwaga: Te 6 par wykazują korelację 1.0, co sugeruje, że są ze sobą liniowo zależne. W modelowaniu ML należy rozważyć usunięcie redundantnych zmiennych, aby uniknąć multikolinearności.
Korelacje bardzo silne (r > 0.70)¶
LUMO ↔ gap → r = +0.8670
LUMO energy jest głównym składnikiem HOMO-LUMO gapelectronic_spatial_extent ↔ U0 → r = +0.8288
Większe cząsteczki mają wyższą energię wewnętrznąpolarizability ↔ U0 → r = +0.7982
Polaryzowalność rośnie z rozmiarem molekułyzpve ↔ U0 → r = +0.7758
Energia drgań punktu zerowego koreluje z całkowitą energiąpolarizability ↔ zpve → r = +0.7655
Obie właściwości zależą od rozmiaru molekułypolarizability ↔ electronic_spatial_extent → r = +0.7170
Bezpośrednia zależność od objętości elektronowej
Korelacje silne (0.50 < r < 0.70)¶
- LUMO ↔ zpve → r = +0.6705
- gap ↔ zpve → r = +0.5705
- electronic_spatial_extent ↔ zpve → r = +0.5284
💡 Kluczowe wnioski:¶
- Właściwości termodynamiczne (U, H, G, Cv) są ze sobą idealnie skorelowane - wystarczy jedna z nich w modelach predykcyjnych
- LUMO i gap są silnie powiązane - gap częściowo zależy od LUMO
- Rozmiar molekuły (mierzony przez polarizability i electronic_spatial_extent) silnie wpływa na właściwości energetyczne
- Brak silnych korelacji ujemnych w top 15 - większość zależności jest dodatnia
Walidacja SMILES - ile błędnych struktur, dlaczego¶
# ============================================================
# KOMÓRKA 7: WALIDACJA I ANALIZA STRUKTUR SMILES
# ============================================================
from rdkit import Chem, RDLogger
from rdkit.Chem import Descriptors, rdMolDescriptors
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
from collections import Counter
# Wyciszenie ostrzeżeń RDKit
RDLogger.DisableLog('rdApp.*')
print("="*70)
print("WALIDACJA STRUKTUR SMILES - QM9 DATASET")
print("="*70)
# ============================================================
# 1. PODSTAWOWA WALIDACJA SMILES
# ============================================================
print("\n📋 ETAP 1: PODSTAWOWA WALIDACJA")
print("-" * 70)
def validate_smiles_detailed(smiles):
"""
Szczegółowa walidacja SMILES z diagnozą błędów
Returns:
dict: {'valid': bool, 'error': str, 'mol': Mol object}
"""
try:
mol = Chem.MolFromSmiles(smiles)
if mol is None:
return {'valid': False, 'error': 'Parse Error', 'mol': None}
# Próba sanityzacji
Chem.SanitizeMol(mol)
return {'valid': True, 'error': None, 'mol': mol}
except Exception as e:
error_type = str(e).split(':')[0] if ':' in str(e) else str(e)
return {'valid': False, 'error': error_type, 'mol': None}
# Walidacja wszystkich SMILES
print("Walidacja struktur SMILES...")
validation_results = df['SMILES'].apply(validate_smiles_detailed)
df['is_valid'] = validation_results.apply(lambda x: x['valid'])
df['error_type'] = validation_results.apply(lambda x: x['error'])
# Statystyki
total_smiles = len(df)
valid_smiles = df['is_valid'].sum()
invalid_smiles = total_smiles - valid_smiles
validity_pct = (valid_smiles / total_smiles) * 100
print(f"\n✅ Poprawne SMILES: {valid_smiles:,} ({validity_pct:.2f}%)")
print(f"❌ Niepoprawne SMILES: {invalid_smiles:,} ({(100-validity_pct):.2f}%)")
# ============================================================
# 2. ANALIZA BŁĘDÓW WALIDACJI
# ============================================================
if invalid_smiles > 0:
print("\n" + "="*70)
print("📊 ANALIZA BŁĘDÓW WALIDACJI")
print("="*70)
# Typy błędów
error_counts = df[~df['is_valid']]['error_type'].value_counts()
print("\n🔍 Typy błędów walidacji:")
for error, count in error_counts.items():
pct = (count / invalid_smiles) * 100
print(f" • {error:40s} {count:5d} ({pct:5.2f}%)")
# Przykłady błędnych SMILES
print("\n⚠️ Przykłady niepoprawnych SMILES:")
invalid_samples = df[~df['is_valid']].head(10)
for idx, row in invalid_samples.iterrows():
print(f"\n [{idx}] {row['SMILES'][:60]}{'...' if len(row['SMILES']) > 60 else ''}")
print(f" Błąd: {row['error_type']}")
# Wizualizacja typów błędów
if len(error_counts) > 0:
plt.figure(figsize=(12, 6))
colors = plt.cm.Reds(np.linspace(0.4, 0.8, len(error_counts)))
plt.barh(range(len(error_counts)), error_counts.values, color=colors, edgecolor='black')
plt.yticks(range(len(error_counts)), error_counts.index)
plt.xlabel('Liczba wystąpień', fontsize=12, fontweight='bold')
plt.title('Typy błędów walidacji SMILES', fontsize=14, fontweight='bold')
plt.grid(axis='x', alpha=0.3)
for i, val in enumerate(error_counts.values):
plt.text(val + max(error_counts.values)*0.01, i, f'{val}',
va='center', fontweight='bold')
plt.tight_layout()
plt.savefig('qm9_smiles_errors.png', dpi=300, bbox_inches='tight')
print("\n✅ Zapisano: qm9_smiles_errors.png")
plt.show()
else:
print("\n✅ WSZYSTKIE SMILES SĄ POPRAWNE! Dataset QM9 ma 100% validację.")
# ============================================================
# 3. ANALIZA SKŁADU CHEMICZNEGO
# ============================================================
print("\n" + "="*70)
print("🧪 ANALIZA SKŁADU CHEMICZNEGO")
print("="*70)
# Tylko dla poprawnych SMILES
df_valid = df[df['is_valid']].copy()
def analyze_composition(smiles):
"""Analiza składu atomowego cząsteczki"""
mol = Chem.MolFromSmiles(smiles)
if mol is None:
return None
atoms = [atom.GetSymbol() for atom in mol.GetAtoms()]
return Counter(atoms)
print("\nAnalizuję skład atomowy cząsteczek...")
compositions = df_valid['SMILES'].apply(analyze_composition)
# Ekstrakcja wszystkich typów atomów
all_atoms = []
for comp in compositions.dropna():
all_atoms.extend(comp.keys())
atom_types = Counter(all_atoms)
print(f"\n📊 Typy atomów występujące w datasecie ({len(atom_types)}):")
for atom, count in atom_types.most_common():
pct = (count / len(df_valid)) * 100
print(f" • {atom:3s} występuje w {count:6,} cząsteczkach ({pct:5.2f}%)")
# Rozkład liczby atomów per typ
atom_counts = {atom: [] for atom in atom_types.keys()}
for comp in compositions.dropna():
for atom in atom_types.keys():
atom_counts[atom].append(comp.get(atom, 0))
# Konwersja do DataFrame
df_atoms = pd.DataFrame(atom_counts)
print(f"\n📈 Statystyki liczby atomów:")
print(df_atoms.describe())
# ============================================================
# 4. WIZUALIZACJA SKŁADU ATOMOWEGO
# ============================================================
print("\n" + "="*70)
print("📊 WIZUALIZACJA SKŁADU ATOMOWEGO")
print("="*70)
# Wykresy dla głównych atomów (H, C, N, O, F)
main_atoms = ['H', 'C', 'N', 'O', 'F']
main_atoms = [atom for atom in main_atoms if atom in df_atoms.columns]
fig, axes = plt.subplots(2, 3, figsize=(18, 10))
axes = axes.flatten()
for i, atom in enumerate(main_atoms):
if atom in df_atoms.columns:
data = df_atoms[atom]
axes[i].hist(data, bins=50, color='steelblue', edgecolor='black', alpha=0.7)
axes[i].set_xlabel(f'Liczba atomów {atom}', fontsize=11, fontweight='bold')
axes[i].set_ylabel('Liczba cząsteczek', fontsize=11, fontweight='bold')
axes[i].set_title(f'Rozkład atomów {atom}\nŚrednia: {data.mean():.1f}, Max: {data.max():.0f}',
fontsize=12, fontweight='bold')
axes[i].grid(axis='y', alpha=0.3)
# Wykres porównawczy - średnie liczby atomów
if len(main_atoms) < 6:
ax_comparison = axes[5]
means = [df_atoms[atom].mean() for atom in main_atoms]
colors = ['#1f77b4', '#ff7f0e', '#2ca02c', '#d62728', '#9467bd'][:len(main_atoms)]
bars = ax_comparison.bar(main_atoms, means, color=colors, edgecolor='black', linewidth=1.5)
ax_comparison.set_ylabel('Średnia liczba atomów', fontsize=11, fontweight='bold')
ax_comparison.set_title('Średni skład atomowy cząsteczek QM9', fontsize=12, fontweight='bold')
ax_comparison.grid(axis='y', alpha=0.3)
for bar, val in zip(bars, means):
ax_comparison.text(bar.get_x() + bar.get_width()/2, val + 0.1,
f'{val:.1f}', ha='center', fontweight='bold')
plt.tight_layout()
plt.savefig('qm9_atomic_composition.png', dpi=300, bbox_inches='tight')
print("✅ Zapisano: qm9_atomic_composition.png")
plt.show()
# ============================================================
# 5. ANALIZA WŁAŚCIWOŚCI MOLEKULARNYCH (RDKit)
# ============================================================
print("\n" + "="*70)
print("⚗️ OBLICZANIE DESKRYPTORÓW MOLEKULARNYCH")
print("="*70)
def calculate_rdkit_descriptors(smiles):
"""Oblicza podstawowe deskryptory RDKit"""
mol = Chem.MolFromSmiles(smiles)
if mol is None:
return None
try:
return {
'MolWt': Descriptors.MolWt(mol),
'NumAtoms': mol.GetNumAtoms(),
'NumHeavyAtoms': mol.GetNumHeavyAtoms(),
'NumHeteroatoms': rdMolDescriptors.CalcNumHeteroatoms(mol),
'NumRotatableBonds': rdMolDescriptors.CalcNumRotatableBonds(mol),
'NumHBD': rdMolDescriptors.CalcNumHBD(mol),
'NumHBA': rdMolDescriptors.CalcNumHBA(mol),
'NumRings': rdMolDescriptors.CalcNumRings(mol),
'NumAromaticRings': rdMolDescriptors.CalcNumAromaticRings(mol),
'TPSA': Descriptors.TPSA(mol),
'LogP': Descriptors.MolLogP(mol),
'FractionCSP3': rdMolDescriptors.CalcFractionCsp3(mol)
}
except:
return None
print("Obliczam deskryptory RDKit dla poprawnych struktur...")
descriptors = df_valid['SMILES'].apply(calculate_rdkit_descriptors)
df_descriptors = descriptors.apply(pd.Series)
# Sprawdzenie czy są jakiekolwiek dane
if df_descriptors.empty or df_descriptors.shape[1] == 0:
print("\n⚠️ BRAK POPRAWNYCH STRUKTUR DO ANALIZY!")
print(" Wszystkie SMILES są niepoprawne lub nie można obliczyć deskryptorów.")
else:
# Dodanie do głównego DataFrame
df_valid = pd.concat([df_valid.reset_index(drop=True),
df_descriptors.reset_index(drop=True)], axis=1)
print(f"✅ Obliczono {len(df_descriptors.columns)} deskryptorów dla {len(df_valid)} struktur")
print("\n📊 Statystyki deskryptorów molekularnych:")
print(df_descriptors.describe())
# ============================================================
# 6. CHARAKTERYSTYKA DATASETU QM9
# ============================================================
if not df_descriptors.empty and df_descriptors.shape[1] > 0:
print("\n" + "="*70)
print("📋 CHARAKTERYSTYKA DATASETU QM9")
print("="*70)
print(f"""
🔬 SKŁAD CHEMICZNY:
• Dozwolone atomy: {', '.join(sorted(atom_types.keys()))}
• Najczęstsze atomy: H, C, N, O, F (QM9 spec)
• Średnia liczba atomów: {df_descriptors['NumAtoms'].mean():.1f}
• Średnia liczba ciężkich atomów: {df_descriptors['NumHeavyAtoms'].mean():.1f}
• Max ciężkich atomów: {df_descriptors['NumHeavyAtoms'].max():.0f} (limit QM9: 9)
📏 ROZMIAR MOLEKUŁ:
• Masa molowa: {df_descriptors['MolWt'].mean():.1f} ± {df_descriptors['MolWt'].std():.1f} g/mol
• Zakres: {df_descriptors['MolWt'].min():.1f} - {df_descriptors['MolWt'].max():.1f} g/mol
🔗 STRUKTURY:
• Średnia liczba pierścieni: {df_descriptors['NumRings'].mean():.2f}
• Cząsteczki aromatyczne: {(df_descriptors['NumAromaticRings'] > 0).sum():,} ({(df_descriptors['NumAromaticRings'] > 0).sum()/len(df_descriptors)*100:.1f}%)
• Średnia rotowalne wiązania: {df_descriptors['NumRotatableBonds'].mean():.2f}
💧 WŁAŚCIWOŚCI FIZYKOCHEMICZNE:
• Średnie LogP: {df_descriptors['LogP'].mean():.2f} ± {df_descriptors['LogP'].std():.2f}
• Średnie TPSA: {df_descriptors['TPSA'].mean():.1f} Ų
• Średnia frakcja sp3: {df_descriptors['FractionCSP3'].mean():.2f}
✅ JAKOŚĆ DANYCH:
• Validacja SMILES: {validity_pct:.2f}%
• Kompletność: {(df['is_valid'].sum() / len(df) * 100):.2f}%
• Błędne struktury: {invalid_smiles:,}
""")
# ============================================================
# 7. WIZUALIZACJA PORÓWNAWCZA
# ============================================================
if not df_descriptors.empty and df_descriptors.shape[1] > 0:
print("\n" + "="*70)
print("📊 WIZUALIZACJA WŁAŚCIWOŚCI MOLEKULARNYCH")
print("="*70)
fig, axes = plt.subplots(2, 3, figsize=(18, 10))
# 1. Masa molowa
axes[0, 0].hist(df_descriptors['MolWt'], bins=50, color='#3498db',
edgecolor='black', alpha=0.7)
axes[0, 0].set_xlabel('Masa molowa (g/mol)', fontweight='bold')
axes[0, 0].set_ylabel('Liczba cząsteczek', fontweight='bold')
axes[0, 0].set_title('Rozkład masy molowej', fontweight='bold', fontsize=12)
axes[0, 0].grid(axis='y', alpha=0.3)
# 2. Liczba ciężkich atomów
axes[0, 1].hist(df_descriptors['NumHeavyAtoms'], bins=range(1, 11),
color='#e74c3c', edgecolor='black', alpha=0.7)
axes[0, 1].set_xlabel('Liczba ciężkich atomów', fontweight='bold')
axes[0, 1].set_ylabel('Liczba cząsteczek', fontweight='bold')
axes[0, 1].set_title('Rozkład liczby ciężkich atomów (max=9)', fontweight='bold', fontsize=12)
axes[0, 1].grid(axis='y', alpha=0.3)
# 3. LogP
axes[0, 2].hist(df_descriptors['LogP'], bins=50, color='#2ecc71',
edgecolor='black', alpha=0.7)
axes[0, 2].set_xlabel('LogP', fontweight='bold')
axes[0, 2].set_ylabel('Liczba cząsteczek', fontweight='bold')
axes[0, 2].set_title('Rozkład lipofilowości (LogP)', fontweight='bold', fontsize=12)
axes[0, 2].grid(axis='y', alpha=0.3)
# 4. TPSA
axes[1, 0].hist(df_descriptors['TPSA'], bins=50, color='#f39c12',
edgecolor='black', alpha=0.7)
axes[1, 0].set_xlabel('TPSA (Ų)', fontweight='bold')
axes[1, 0].set_ylabel('Liczba cząsteczek', fontweight='bold')
axes[1, 0].set_title('Rozkład polarnej powierzchni (TPSA)', fontweight='bold', fontsize=12)
axes[1, 0].grid(axis='y', alpha=0.3)
# 5. Liczba pierścieni
ring_counts = df_descriptors['NumRings'].value_counts().sort_index()
axes[1, 1].bar(ring_counts.index, ring_counts.values, color='#9b59b6',
edgecolor='black', alpha=0.7)
axes[1, 1].set_xlabel('Liczba pierścieni', fontweight='bold')
axes[1, 1].set_ylabel('Liczba cząsteczek', fontweight='bold')
axes[1, 1].set_title('Rozkład liczby pierścieni', fontweight='bold', fontsize=12)
axes[1, 1].grid(axis='y', alpha=0.3)
# 6. Frakcja sp3
axes[1, 2].hist(df_descriptors['FractionCSP3'], bins=50, color='#1abc9c',
edgecolor='black', alpha=0.7)
axes[1, 2].set_xlabel('Frakcja węgli sp³', fontweight='bold')
axes[1, 2].set_ylabel('Liczba cząsteczek', fontweight='bold')
axes[1, 2].set_title('Rozkład frakcji sp³ (saturacja)', fontweight='bold', fontsize=12)
axes[1, 2].grid(axis='y', alpha=0.3)
plt.suptitle('Właściwości molekularne datasetu QM9', fontsize=16, fontweight='bold', y=1.00)
plt.tight_layout()
plt.savefig('qm9_molecular_properties.png', dpi=300, bbox_inches='tight')
print("✅ Zapisano: qm9_molecular_properties.png")
plt.show()
else:
print("\n⚠️ Brak danych do wizualizacji właściwości molekularnych")
print("\n" + "="*70)
print("✅ WALIDACJA I ANALIZA SMILES ZAKOŃCZONA!")
print("="*70)
======================================================================
WALIDACJA STRUKTUR SMILES - QM9 DATASET
======================================================================
📋 ETAP 1: PODSTAWOWA WALIDACJA
----------------------------------------------------------------------
Walidacja struktur SMILES...
✅ Poprawne SMILES: 1,000 (100.00%)
❌ Niepoprawne SMILES: 0 (0.00%)
✅ WSZYSTKIE SMILES SĄ POPRAWNE! Dataset QM9 ma 100% validację.
======================================================================
🧪 ANALIZA SKŁADU CHEMICZNEGO
======================================================================
Analizuję skład atomowy cząsteczek...
📊 Typy atomów występujące w datasecie (5):
• C występuje w 1,000 cząsteczkach (100.00%)
• O występuje w 843 cząsteczkach (84.30%)
• N występuje w 622 cząsteczkach (62.20%)
• H występuje w 16 cząsteczkach ( 1.60%)
• F występuje w 14 cząsteczkach ( 1.40%)
📈 Statystyki liczby atomów:
C N O H F
count 1000.000000 1000.000000 1000.000000 1000.000000 1000.000000
mean 6.251000 1.093000 1.392000 0.017000 0.020000
std 1.282046 1.138265 0.875844 0.136857 0.188774
min 2.000000 0.000000 0.000000 0.000000 0.000000
25% 5.000000 0.000000 1.000000 0.000000 0.000000
50% 6.000000 1.000000 1.000000 0.000000 0.000000
75% 7.000000 2.000000 2.000000 0.000000 0.000000
max 9.000000 6.000000 4.000000 2.000000 3.000000
======================================================================
📊 WIZUALIZACJA SKŁADU ATOMOWEGO
======================================================================
✅ Zapisano: qm9_atomic_composition.png
====================================================================== ⚗️ OBLICZANIE DESKRYPTORÓW MOLEKULARNYCH ====================================================================== Obliczam deskryptory RDKit dla poprawnych struktur... ⚠️ BRAK POPRAWNYCH STRUKTUR DO ANALIZY! Wszystkie SMILES są niepoprawne lub nie można obliczyć deskryptorów. ⚠️ Brak danych do wizualizacji właściwości molekularnych ====================================================================== ✅ WALIDACJA I ANALIZA SMILES ZAKOŃCZONA! ======================================================================
# debudding
# ============================================================
# DEBUGGING - SPRAWDZENIE KOLUMNY SMILES
# ============================================================
print("="*70)
print("🔍 DEBUGGING - ANALIZA KOLUMNY SMILES")
print("="*70)
# 1. Sprawdź nazwy kolumn
print("\n1️⃣ Kolumny w DataFrame:")
print(df.columns.tolist())
# 2. Sprawdź czy kolumna SMILES istnieje
print(f"\n2️⃣ Czy 'SMILES' w kolumnach? {('SMILES' in df.columns)}")
# 3. Jeśli istnieje, pokaż przykłady
if 'SMILES' in df.columns:
print("\n3️⃣ Pierwsze 10 wartości w kolumnie SMILES:")
for i, smiles in enumerate(df['SMILES'].head(10), 1):
print(f" {i:2d}. {smiles}")
# 4. Typ danych
print(f"\n4️⃣ Typ danych kolumny SMILES: {df['SMILES'].dtype}")
# 5. Długość SMILES
print(f"\n5️⃣ Statystyki długości SMILES:")
smiles_lengths = df['SMILES'].astype(str).str.len()
print(f" Min: {smiles_lengths.min()}")
print(f" Max: {smiles_lengths.max()}")
print(f" Średnia: {smiles_lengths.mean():.1f}")
# 6. Czy są puste wartości?
print(f"\n6️⃣ Puste/NaN wartości: {df['SMILES'].isna().sum()}")
# 7. Test walidacji pojedynczego SMILES
print("\n7️⃣ Test walidacji pojedynczego SMILES:")
test_smiles = df['SMILES'].iloc[0]
print(f" SMILES: {test_smiles}")
print(f" Typ: {type(test_smiles)}")
from rdkit import Chem
mol = Chem.MolFromSmiles(str(test_smiles))
if mol is None:
print(f" ❌ RDKit NIE MOŻE SPARSOWAĆ!")
else:
print(f" ✅ RDKit sparsował poprawnie")
print(f" Liczba atomów: {mol.GetNumAtoms()}")
else:
print("\n❌ KOLUMNA 'SMILES' NIE ISTNIEJE!")
print("\nMożliwe alternatywne nazwy:")
for col in df.columns:
if 'smile' in col.lower() or 'mol' in col.lower() or 'id' in col.lower():
print(f" • {col}")
print("\n" + "="*70)
====================================================================== 🔍 DEBUGGING - ANALIZA KOLUMNY SMILES ====================================================================== 1️⃣ Kolumny w DataFrame: ['dipole_moment', 'polarizability', 'HOMO', 'LUMO', 'gap', 'electronic_spatial_extent', 'zpve', 'U0', 'U', 'H', 'G', 'Cv', 'SMILES', 'dataset_split', 'num_atoms', 'is_valid', 'error_type'] 2️⃣ Czy 'SMILES' w kolumnach? True 3️⃣ Pierwsze 10 wartości w kolumnie SMILES: 1. [H]C([H])([H])C([H])([H])[N@@H+]1C([H])([H])[C@]1(C([H])([H])[H])C([H])([H])C([H])([H])C([H])([H])[H] 2. [H]O[C@@]12C([H])([H])O[C@]1([H])[C@]2([H])[C@]([H])(O[H])C([H])([H])[H] 3. [H]C1([H])C([H])([H])[C@]2([H])[C@@]1([H])[C@]1([H])OC([H])([H])[C@@]12C([H])([H])[H] 4. [H]OC([H])([H])[C@]1(C([H])([H])C#N)OC1([H])[H] 5. [H]N1[C@]2([H])C(=O)[C@@]1([H])C([H])([H])OC2([H])[H] 6. [H][N-]C1OC([H])([H])[C@]12[N@@H+]1C([H])([H])[C@]1([H])C2([H])[H] 7. [H]N([H])C1[NH2+][C@]2([H])C([H])([H])[C@]2(C([H])([H])C([H])([H])[H])O1 8. [H]C1([H])[N@@H+]2[C@@]3([H])[C@]2([H])[C@@]2([H])C([H])([H])[C@]1([H])[C@@]32[H] 9. [H]C([H])([H])O[C@@]1(C([H])([H])[H])C([H])([H])[C@@]2([H])[C@@]([H])(C([H])([H])[H])[C@@]12[H] 10. [H]C(=O)c1c([H])c(N([H])[H])nn1[H] 4️⃣ Typ danych kolumny SMILES: object 5️⃣ Statystyki długości SMILES: Min: 13 Max: 101 Średnia: 62.8 6️⃣ Puste/NaN wartości: 0 7️⃣ Test walidacji pojedynczego SMILES: SMILES: [H]C([H])([H])C([H])([H])[N@@H+]1C([H])([H])[C@]1(C([H])([H])[H])C([H])([H])C([H])([H])C([H])([H])[H] Typ: <class 'str'> ✅ RDKit sparsował poprawnie Liczba atomów: 9 ======================================================================
SMILES są poprawne, ale mają jawne atomy wodoru [H] i sterochemię [C@], [C@@] - to są bardzo szczegółowe SMILES z QM9! Problem jest w poprzednim kodzie - prawdopodobnie kolumny is_valid i error_type już istnieją z poprzedniego uruchomienia i nadpisują nowe wyniki. Rozwiązanie: Usuńmy te kolumny przed walidacją:
# ============================================================
# KOMÓRKA 7: WALIDACJA I ANALIZA STRUKTUR SMILES
# ============================================================
from rdkit import Chem, RDLogger
from rdkit.Chem import Descriptors, rdMolDescriptors
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
from collections import Counter
# Wyciszenie ostrzeżeń RDKit
RDLogger.DisableLog('rdApp.*')
print("="*70)
print("WALIDACJA STRUKTUR SMILES - QM9 DATASET")
print("="*70)
# ============================================================
# 1. PODSTAWOWA WALIDACJA SMILES
# ============================================================
print("\n📋 ETAP 1: PODSTAWOWA WALIDACJA")
print("-" * 70)
# Usunięcie starych kolumn walidacji jeśli istnieją
if 'is_valid' in df.columns:
df = df.drop(columns=['is_valid'])
if 'error_type' in df.columns:
df = df.drop(columns=['error_type'])
def validate_smiles_detailed(smiles):
"""
Szczegółowa walidacja SMILES z diagnozą błędów
Returns:
dict: {'valid': bool, 'error': str, 'mol': Mol object}
"""
try:
mol = Chem.MolFromSmiles(smiles)
if mol is None:
return {'valid': False, 'error': 'Parse Error', 'mol': None}
# Próba sanityzacji
Chem.SanitizeMol(mol)
return {'valid': True, 'error': None, 'mol': mol}
except Exception as e:
error_type = str(e).split(':')[0] if ':' in str(e) else str(e)
return {'valid': False, 'error': error_type, 'mol': None}
# Walidacja wszystkich SMILES
print("Walidacja struktur SMILES...")
print(f"Przykładowe SMILES do walidacji:")
print(df['SMILES'].head(5).tolist())
validation_results = df['SMILES'].apply(validate_smiles_detailed)
df['is_valid'] = validation_results.apply(lambda x: x['valid'])
df['error_type'] = validation_results.apply(lambda x: x['error'])
# Statystyki
total_smiles = len(df)
valid_smiles = df['is_valid'].sum()
invalid_smiles = total_smiles - valid_smiles
validity_pct = (valid_smiles / total_smiles) * 100
print(f"\n✅ Poprawne SMILES: {valid_smiles:,} ({validity_pct:.2f}%)")
print(f"❌ Niepoprawne SMILES: {invalid_smiles:,} ({(100-validity_pct):.2f}%)")
# Debug - pokaż przykłady
if valid_smiles > 0:
print(f"\n✓ Przykłady POPRAWNYCH SMILES:")
print(df[df['is_valid']]['SMILES'].head(3).tolist())
if invalid_smiles > 0:
print(f"\n✗ Przykłady NIEPOPRAWNYCH SMILES:")
print(df[~df['is_valid']]['SMILES'].head(3).tolist())
# ============================================================
# 2. ANALIZA BŁĘDÓW WALIDACJI
# ============================================================
if invalid_smiles > 0:
print("\n" + "="*70)
print("📊 ANALIZA BŁĘDÓW WALIDACJI")
print("="*70)
# Typy błędów
error_counts = df[~df['is_valid']]['error_type'].value_counts()
print("\n🔍 Typy błędów walidacji:")
for error, count in error_counts.items():
pct = (count / invalid_smiles) * 100
print(f" • {error:40s} {count:5d} ({pct:5.2f}%)")
# Przykłady błędnych SMILES
print("\n⚠️ Przykłady niepoprawnych SMILES:")
invalid_samples = df[~df['is_valid']].head(10)
for idx, row in invalid_samples.iterrows():
print(f"\n [{idx}] {row['SMILES'][:60]}{'...' if len(row['SMILES']) > 60 else ''}")
print(f" Błąd: {row['error_type']}")
# Wizualizacja typów błędów
if len(error_counts) > 0:
plt.figure(figsize=(12, 6))
colors = plt.cm.Reds(np.linspace(0.4, 0.8, len(error_counts)))
plt.barh(range(len(error_counts)), error_counts.values, color=colors, edgecolor='black')
plt.yticks(range(len(error_counts)), error_counts.index)
plt.xlabel('Liczba wystąpień', fontsize=12, fontweight='bold')
plt.title('Typy błędów walidacji SMILES', fontsize=14, fontweight='bold')
plt.grid(axis='x', alpha=0.3)
for i, val in enumerate(error_counts.values):
plt.text(val + max(error_counts.values)*0.01, i, f'{val}',
va='center', fontweight='bold')
plt.tight_layout()
plt.savefig('qm9_smiles_errors.png', dpi=300, bbox_inches='tight')
print("\n✅ Zapisano: qm9_smiles_errors.png")
plt.show()
else:
print("\n✅ WSZYSTKIE SMILES SĄ POPRAWNE! Dataset QM9 ma 100% validację.")
# ============================================================
# 3. ANALIZA SKŁADU CHEMICZNEGO
# ============================================================
print("\n" + "="*70)
print("🧪 ANALIZA SKŁADU CHEMICZNEGO")
print("="*70)
# Tylko dla poprawnych SMILES
df_valid = df[df['is_valid']].copy()
def analyze_composition(smiles):
"""Analiza składu atomowego cząsteczki"""
mol = Chem.MolFromSmiles(smiles)
if mol is None:
return None
atoms = [atom.GetSymbol() for atom in mol.GetAtoms()]
return Counter(atoms)
print("\nAnalizuję skład atomowy cząsteczek...")
compositions = df_valid['SMILES'].apply(analyze_composition)
# Ekstrakcja wszystkich typów atomów
all_atoms = []
for comp in compositions.dropna():
all_atoms.extend(comp.keys())
atom_types = Counter(all_atoms)
print(f"\n📊 Typy atomów występujące w datasecie ({len(atom_types)}):")
for atom, count in atom_types.most_common():
pct = (count / len(df_valid)) * 100
print(f" • {atom:3s} występuje w {count:6,} cząsteczkach ({pct:5.2f}%)")
# Rozkład liczby atomów per typ
atom_counts = {atom: [] for atom in atom_types.keys()}
for comp in compositions.dropna():
for atom in atom_types.keys():
atom_counts[atom].append(comp.get(atom, 0))
# Konwersja do DataFrame
df_atoms = pd.DataFrame(atom_counts)
print(f"\n📈 Statystyki liczby atomów:")
print(df_atoms.describe())
# ============================================================
# 4. WIZUALIZACJA SKŁADU ATOMOWEGO
# ============================================================
print("\n" + "="*70)
print("📊 WIZUALIZACJA SKŁADU ATOMOWEGO")
print("="*70)
# Wykresy dla głównych atomów (H, C, N, O, F)
main_atoms = ['H', 'C', 'N', 'O', 'F']
main_atoms = [atom for atom in main_atoms if atom in df_atoms.columns]
fig, axes = plt.subplots(2, 3, figsize=(18, 10))
axes = axes.flatten()
for i, atom in enumerate(main_atoms):
if atom in df_atoms.columns:
data = df_atoms[atom]
axes[i].hist(data, bins=50, color='steelblue', edgecolor='black', alpha=0.7)
axes[i].set_xlabel(f'Liczba atomów {atom}', fontsize=11, fontweight='bold')
axes[i].set_ylabel('Liczba cząsteczek', fontsize=11, fontweight='bold')
axes[i].set_title(f'Rozkład atomów {atom}\nŚrednia: {data.mean():.1f}, Max: {data.max():.0f}',
fontsize=12, fontweight='bold')
axes[i].grid(axis='y', alpha=0.3)
# Wykres porównawczy - średnie liczby atomów
if len(main_atoms) < 6:
ax_comparison = axes[5]
means = [df_atoms[atom].mean() for atom in main_atoms]
colors = ['#1f77b4', '#ff7f0e', '#2ca02c', '#d62728', '#9467bd'][:len(main_atoms)]
bars = ax_comparison.bar(main_atoms, means, color=colors, edgecolor='black', linewidth=1.5)
ax_comparison.set_ylabel('Średnia liczba atomów', fontsize=11, fontweight='bold')
ax_comparison.set_title('Średni skład atomowy cząsteczek QM9', fontsize=12, fontweight='bold')
ax_comparison.grid(axis='y', alpha=0.3)
for bar, val in zip(bars, means):
ax_comparison.text(bar.get_x() + bar.get_width()/2, val + 0.1,
f'{val:.1f}', ha='center', fontweight='bold')
plt.tight_layout()
plt.savefig('qm9_atomic_composition.png', dpi=300, bbox_inches='tight')
print("✅ Zapisano: qm9_atomic_composition.png")
plt.show()
# ============================================================
# 5. ANALIZA WŁAŚCIWOŚCI MOLEKULARNYCH (RDKit)
# ============================================================
print("\n" + "="*70)
print("⚗️ OBLICZANIE DESKRYPTORÓW MOLEKULARNYCH")
print("="*70)
def calculate_rdkit_descriptors(smiles):
"""Oblicza podstawowe deskryptory RDKit"""
mol = Chem.MolFromSmiles(smiles)
if mol is None:
return None
try:
return {
'MolWt': Descriptors.MolWt(mol),
'NumAtoms': mol.GetNumAtoms(),
'NumHeavyAtoms': mol.GetNumHeavyAtoms(),
'NumHeteroatoms': rdMolDescriptors.CalcNumHeteroatoms(mol),
'NumRotatableBonds': rdMolDescriptors.CalcNumRotatableBonds(mol),
'NumHBD': rdMolDescriptors.CalcNumHBD(mol),
'NumHBA': rdMolDescriptors.CalcNumHBA(mol),
'NumRings': rdMolDescriptors.CalcNumRings(mol),
'NumAromaticRings': rdMolDescriptors.CalcNumAromaticRings(mol),
'TPSA': Descriptors.TPSA(mol),
'LogP': Descriptors.MolLogP(mol),
'FractionCSP3': rdMolDescriptors.CalcFractionCsp3(mol)
}
except:
return None
print("Obliczam deskryptory RDKit dla poprawnych struktur...")
descriptors = df_valid['SMILES'].apply(calculate_rdkit_descriptors)
df_descriptors = descriptors.apply(pd.Series)
# Sprawdzenie czy są jakiekolwiek dane
if df_descriptors.empty or df_descriptors.shape[1] == 0:
print("\n⚠️ BRAK POPRAWNYCH STRUKTUR DO ANALIZY!")
print(" Wszystkie SMILES są niepoprawne lub nie można obliczyć deskryptorów.")
else:
# Dodanie do głównego DataFrame
df_valid = pd.concat([df_valid.reset_index(drop=True),
df_descriptors.reset_index(drop=True)], axis=1)
print(f"✅ Obliczono {len(df_descriptors.columns)} deskryptorów dla {len(df_valid)} struktur")
print("\n📊 Statystyki deskryptorów molekularnych:")
print(df_descriptors.describe())
# ============================================================
# 6. CHARAKTERYSTYKA DATASETU QM9
# ============================================================
if not df_descriptors.empty and df_descriptors.shape[1] > 0:
print("\n" + "="*70)
print("📋 CHARAKTERYSTYKA DATASETU QM9")
print("="*70)
print(f"""
🔬 SKŁAD CHEMICZNY:
• Dozwolone atomy: {', '.join(sorted(atom_types.keys()))}
• Najczęstsze atomy: H, C, N, O, F (QM9 spec)
• Średnia liczba atomów: {df_descriptors['NumAtoms'].mean():.1f}
• Średnia liczba ciężkich atomów: {df_descriptors['NumHeavyAtoms'].mean():.1f}
• Max ciężkich atomów: {df_descriptors['NumHeavyAtoms'].max():.0f} (limit QM9: 9)
📏 ROZMIAR MOLEKUŁ:
• Masa molowa: {df_descriptors['MolWt'].mean():.1f} ± {df_descriptors['MolWt'].std():.1f} g/mol
• Zakres: {df_descriptors['MolWt'].min():.1f} - {df_descriptors['MolWt'].max():.1f} g/mol
🔗 STRUKTURY:
• Średnia liczba pierścieni: {df_descriptors['NumRings'].mean():.2f}
• Cząsteczki aromatyczne: {(df_descriptors['NumAromaticRings'] > 0).sum():,} ({(df_descriptors['NumAromaticRings'] > 0).sum()/len(df_descriptors)*100:.1f}%)
• Średnia rotowalne wiązania: {df_descriptors['NumRotatableBonds'].mean():.2f}
💧 WŁAŚCIWOŚCI FIZYKOCHEMICZNE:
• Średnie LogP: {df_descriptors['LogP'].mean():.2f} ± {df_descriptors['LogP'].std():.2f}
• Średnie TPSA: {df_descriptors['TPSA'].mean():.1f} Ų
• Średnia frakcja sp3: {df_descriptors['FractionCSP3'].mean():.2f}
✅ JAKOŚĆ DANYCH:
• Validacja SMILES: {validity_pct:.2f}%
• Kompletność: {(df['is_valid'].sum() / len(df) * 100):.2f}%
• Błędne struktury: {invalid_smiles:,}
""")
# ============================================================
# 7. WIZUALIZACJA PORÓWNAWCZA
# ============================================================
if not df_descriptors.empty and df_descriptors.shape[1] > 0:
print("\n" + "="*70)
print("📊 WIZUALIZACJA WŁAŚCIWOŚCI MOLEKULARNYCH")
print("="*70)
fig, axes = plt.subplots(2, 3, figsize=(18, 10))
# 1. Masa molowa
axes[0, 0].hist(df_descriptors['MolWt'], bins=50, color='#3498db',
edgecolor='black', alpha=0.7)
axes[0, 0].set_xlabel('Masa molowa (g/mol)', fontweight='bold')
axes[0, 0].set_ylabel('Liczba cząsteczek', fontweight='bold')
axes[0, 0].set_title('Rozkład masy molowej', fontweight='bold', fontsize=12)
axes[0, 0].grid(axis='y', alpha=0.3)
# 2. Liczba ciężkich atomów
axes[0, 1].hist(df_descriptors['NumHeavyAtoms'], bins=range(1, 11),
color='#e74c3c', edgecolor='black', alpha=0.7)
axes[0, 1].set_xlabel('Liczba ciężkich atomów', fontweight='bold')
axes[0, 1].set_ylabel('Liczba cząsteczek', fontweight='bold')
axes[0, 1].set_title('Rozkład liczby ciężkich atomów (max=9)', fontweight='bold', fontsize=12)
axes[0, 1].grid(axis='y', alpha=0.3)
# 3. LogP
axes[0, 2].hist(df_descriptors['LogP'], bins=50, color='#2ecc71',
edgecolor='black', alpha=0.7)
axes[0, 2].set_xlabel('LogP', fontweight='bold')
axes[0, 2].set_ylabel('Liczba cząsteczek', fontweight='bold')
axes[0, 2].set_title('Rozkład lipofilowości (LogP)', fontweight='bold', fontsize=12)
axes[0, 2].grid(axis='y', alpha=0.3)
# 4. TPSA
axes[1, 0].hist(df_descriptors['TPSA'], bins=50, color='#f39c12',
edgecolor='black', alpha=0.7)
axes[1, 0].set_xlabel('TPSA (Ų)', fontweight='bold')
axes[1, 0].set_ylabel('Liczba cząsteczek', fontweight='bold')
axes[1, 0].set_title('Rozkład polarnej powierzchni (TPSA)', fontweight='bold', fontsize=12)
axes[1, 0].grid(axis='y', alpha=0.3)
# 5. Liczba pierścieni
ring_counts = df_descriptors['NumRings'].value_counts().sort_index()
axes[1, 1].bar(ring_counts.index, ring_counts.values, color='#9b59b6',
edgecolor='black', alpha=0.7)
axes[1, 1].set_xlabel('Liczba pierścieni', fontweight='bold')
axes[1, 1].set_ylabel('Liczba cząsteczek', fontweight='bold')
axes[1, 1].set_title('Rozkład liczby pierścieni', fontweight='bold', fontsize=12)
axes[1, 1].grid(axis='y', alpha=0.3)
# 6. Frakcja sp3
axes[1, 2].hist(df_descriptors['FractionCSP3'], bins=50, color='#1abc9c',
edgecolor='black', alpha=0.7)
axes[1, 2].set_xlabel('Frakcja węgli sp³', fontweight='bold')
axes[1, 2].set_ylabel('Liczba cząsteczek', fontweight='bold')
axes[1, 2].set_title('Rozkład frakcji sp³ (saturacja)', fontweight='bold', fontsize=12)
axes[1, 2].grid(axis='y', alpha=0.3)
plt.suptitle('Właściwości molekularne datasetu QM9', fontsize=16, fontweight='bold', y=1.00)
plt.tight_layout()
plt.savefig('qm9_molecular_properties.png', dpi=300, bbox_inches='tight')
print("✅ Zapisano: qm9_molecular_properties.png")
plt.show()
else:
print("\n⚠️ Brak danych do wizualizacji właściwości molekularnych")
print("\n" + "="*70)
print("✅ WALIDACJA I ANALIZA SMILES ZAKOŃCZONA!")
print("="*70)
======================================================================
WALIDACJA STRUKTUR SMILES - QM9 DATASET
======================================================================
📋 ETAP 1: PODSTAWOWA WALIDACJA
----------------------------------------------------------------------
Walidacja struktur SMILES...
Przykładowe SMILES do walidacji:
['[H]C([H])([H])C([H])([H])[N@@H+]1C([H])([H])[C@]1(C([H])([H])[H])C([H])([H])C([H])([H])C([H])([H])[H]', '[H]O[C@@]12C([H])([H])O[C@]1([H])[C@]2([H])[C@]([H])(O[H])C([H])([H])[H]', '[H]C1([H])C([H])([H])[C@]2([H])[C@@]1([H])[C@]1([H])OC([H])([H])[C@@]12C([H])([H])[H]', '[H]OC([H])([H])[C@]1(C([H])([H])C#N)OC1([H])[H]', '[H]N1[C@]2([H])C(=O)[C@@]1([H])C([H])([H])OC2([H])[H]']
✅ Poprawne SMILES: 1,000 (100.00%)
❌ Niepoprawne SMILES: 0 (0.00%)
✓ Przykłady POPRAWNYCH SMILES:
['[H]C([H])([H])C([H])([H])[N@@H+]1C([H])([H])[C@]1(C([H])([H])[H])C([H])([H])C([H])([H])C([H])([H])[H]', '[H]O[C@@]12C([H])([H])O[C@]1([H])[C@]2([H])[C@]([H])(O[H])C([H])([H])[H]', '[H]C1([H])C([H])([H])[C@]2([H])[C@@]1([H])[C@]1([H])OC([H])([H])[C@@]12C([H])([H])[H]']
✅ WSZYSTKIE SMILES SĄ POPRAWNE! Dataset QM9 ma 100% validację.
======================================================================
🧪 ANALIZA SKŁADU CHEMICZNEGO
======================================================================
Analizuję skład atomowy cząsteczek...
📊 Typy atomów występujące w datasecie (5):
• C występuje w 1,000 cząsteczkach (100.00%)
• O występuje w 843 cząsteczkach (84.30%)
• N występuje w 622 cząsteczkach (62.20%)
• H występuje w 16 cząsteczkach ( 1.60%)
• F występuje w 14 cząsteczkach ( 1.40%)
📈 Statystyki liczby atomów:
C N O H F
count 1000.000000 1000.000000 1000.000000 1000.000000 1000.000000
mean 6.251000 1.093000 1.392000 0.017000 0.020000
std 1.282046 1.138265 0.875844 0.136857 0.188774
min 2.000000 0.000000 0.000000 0.000000 0.000000
25% 5.000000 0.000000 1.000000 0.000000 0.000000
50% 6.000000 1.000000 1.000000 0.000000 0.000000
75% 7.000000 2.000000 2.000000 0.000000 0.000000
max 9.000000 6.000000 4.000000 2.000000 3.000000
======================================================================
📊 WIZUALIZACJA SKŁADU ATOMOWEGO
======================================================================
✅ Zapisano: qm9_atomic_composition.png
====================================================================== ⚗️ OBLICZANIE DESKRYPTORÓW MOLEKULARNYCH ====================================================================== Obliczam deskryptory RDKit dla poprawnych struktur... ⚠️ BRAK POPRAWNYCH STRUKTUR DO ANALIZY! Wszystkie SMILES są niepoprawne lub nie można obliczyć deskryptorów. ⚠️ Brak danych do wizualizacji właściwości molekularnych ====================================================================== ✅ WALIDACJA I ANALIZA SMILES ZAKOŃCZONA! ======================================================================
# ============================================================
# KOMÓRKA 7: WALIDACJA I ANALIZA STRUKTUR SMILES
# ============================================================
from rdkit import Chem, RDLogger
from rdkit.Chem import Descriptors, rdMolDescriptors
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
from collections import Counter
# Wyciszenie ostrzeżeń RDKit
RDLogger.DisableLog('rdApp.*')
print("="*70)
print("WALIDACJA STRUKTUR SMILES - QM9 DATASET")
print("="*70)
# ============================================================
# 1. PODSTAWOWA WALIDACJA SMILES
# ============================================================
print("\n📋 ETAP 1: PODSTAWOWA WALIDACJA")
print("-" * 70)
# Usunięcie starych kolumn walidacji jeśli istnieją
if 'is_valid' in df.columns:
df = df.drop(columns=['is_valid'])
if 'error_type' in df.columns:
df = df.drop(columns=['error_type'])
def validate_smiles_detailed(smiles):
"""Szczegółowa walidacja SMILES z diagnozą błędów"""
try:
mol = Chem.MolFromSmiles(smiles)
if mol is None:
return {'valid': False, 'error': 'Parse Error', 'mol': None}
# Próba sanityzacji
Chem.SanitizeMol(mol)
return {'valid': True, 'error': None, 'mol': mol}
except Exception as e:
error_type = str(e).split(':')[0] if ':' in str(e) else str(e)
return {'valid': False, 'error': error_type, 'mol': None}
# Walidacja wszystkich SMILES
print("Walidacja struktur SMILES...")
validation_results = df['SMILES'].apply(validate_smiles_detailed)
df['is_valid'] = validation_results.apply(lambda x: x['valid'])
df['error_type'] = validation_results.apply(lambda x: x['error'])
# Statystyki
total_smiles = len(df)
valid_smiles = df['is_valid'].sum()
invalid_smiles = total_smiles - valid_smiles
validity_pct = (valid_smiles / total_smiles) * 100
print(f"\n✅ Poprawne SMILES: {valid_smiles:,} ({validity_pct:.2f}%)")
print(f"❌ Niepoprawne SMILES: {invalid_smiles:,} ({(100-validity_pct):.2f}%)")
# ============================================================
# 2. SPRAWDZENIE CZY MAMY POPRAWNE DANE
# ============================================================
if valid_smiles == 0:
print("\n" + "="*70)
print("❌ BŁĄD: BRAK POPRAWNYCH STRUKTUR!")
print("="*70)
print("\nPrawdopodobne przyczyny:")
print(" 1. Kolumna 'is_valid' miała już wartości False z poprzedniego uruchomienia")
print(" 2. Problem z biblioteką RDKit")
print("\nSpróbuj zrestartować kernel i uruchomić ponownie od początku.")
else:
print(f"\n✓ Kontynuuję analizę {valid_smiles:,} poprawnych struktur...")
# Tylko dla poprawnych SMILES
df_valid = df[df['is_valid']].copy()
# ============================================================
# 3. ANALIZA SKŁADU CHEMICZNEGO
# ============================================================
if len(df_valid) > 0:
print("\n" + "="*70)
print("🧪 ANALIZA SKŁADU CHEMICZNEGO")
print("="*70)
def analyze_composition(smiles):
"""Analiza składu atomowego cząsteczki"""
mol = Chem.MolFromSmiles(smiles)
if mol is None:
return None
atoms = [atom.GetSymbol() for atom in mol.GetAtoms()]
return Counter(atoms)
print("\nAnalizuję skład atomowy cząsteczek...")
compositions = df_valid['SMILES'].apply(analyze_composition)
# Ekstrakcja wszystkich typów atomów
all_atoms = []
for comp in compositions.dropna():
all_atoms.extend(comp.keys())
atom_types = Counter(all_atoms)
print(f"\n📊 Typy atomów występujące w datasecie ({len(atom_types)}):")
for atom, count in atom_types.most_common():
pct = (count / len(df_valid)) * 100
print(f" • {atom:3s} występuje w {count:6,} cząsteczkach ({pct:5.2f}%)")
# Rozkład liczby atomów per typ
atom_counts = {atom: [] for atom in atom_types.keys()}
for comp in compositions.dropna():
for atom in atom_types.keys():
atom_counts[atom].append(comp.get(atom, 0))
# Konwersja do DataFrame
df_atoms = pd.DataFrame(atom_counts)
print(f"\n📈 Statystyki liczby atomów:")
print(df_atoms.describe())
# ============================================================
# 4. WIZUALIZACJA SKŁADU ATOMOWEGO
# ============================================================
print("\n" + "="*70)
print("📊 WIZUALIZACJA SKŁADU ATOMOWEGO")
print("="*70)
# Wykresy dla głównych atomów (H, C, N, O, F)
main_atoms = ['H', 'C', 'N', 'O', 'F']
main_atoms = [atom for atom in main_atoms if atom in df_atoms.columns]
fig, axes = plt.subplots(2, 3, figsize=(18, 10))
axes = axes.flatten()
for i, atom in enumerate(main_atoms):
if atom in df_atoms.columns:
data = df_atoms[atom]
axes[i].hist(data, bins=50, color='steelblue', edgecolor='black', alpha=0.7)
axes[i].set_xlabel(f'Liczba atomów {atom}', fontsize=11, fontweight='bold')
axes[i].set_ylabel('Liczba cząsteczek', fontsize=11, fontweight='bold')
axes[i].set_title(f'Rozkład atomów {atom}\nŚrednia: {data.mean():.1f}, Max: {data.max():.0f}',
fontsize=12, fontweight='bold')
axes[i].grid(axis='y', alpha=0.3)
# Wykres porównawczy - średnie liczby atomów
if len(main_atoms) < 6:
ax_comparison = axes[5]
means = [df_atoms[atom].mean() for atom in main_atoms]
colors = ['#1f77b4', '#ff7f0e', '#2ca02c', '#d62728', '#9467bd'][:len(main_atoms)]
bars = ax_comparison.bar(main_atoms, means, color=colors, edgecolor='black', linewidth=1.5)
ax_comparison.set_ylabel('Średnia liczba atomów', fontsize=11, fontweight='bold')
ax_comparison.set_title('Średni skład atomowy cząsteczek QM9', fontsize=12, fontweight='bold')
ax_comparison.grid(axis='y', alpha=0.3)
for bar, val in zip(bars, means):
ax_comparison.text(bar.get_x() + bar.get_width()/2, val + 0.1,
f'{val:.1f}', ha='center', fontweight='bold')
plt.tight_layout()
plt.savefig('qm9_atomic_composition.png', dpi=300, bbox_inches='tight')
print("✅ Zapisano: qm9_atomic_composition.png")
plt.show()
# ============================================================
# 5. ANALIZA WŁAŚCIWOŚCI MOLEKULARNYCH (RDKit)
# ============================================================
print("\n" + "="*70)
print("⚗️ OBLICZANIE DESKRYPTORÓW MOLEKULARNYCH")
print("="*70)
def calculate_rdkit_descriptors(smiles):
"""Oblicza podstawowe deskryptory RDKit"""
mol = Chem.MolFromSmiles(smiles)
if mol is None:
return None
try:
return {
'MolWt': Descriptors.MolWt(mol),
'NumAtoms': mol.GetNumAtoms(),
'NumHeavyAtoms': mol.GetNumHeavyAtoms(),
'NumHeteroatoms': rdMolDescriptors.CalcNumHeteroatoms(mol),
'NumRotatableBonds': rdMolDescriptors.CalcNumRotatableBonds(mol),
'NumHBD': rdMolDescriptors.CalcNumHBD(mol),
'NumHBA': rdMolDescriptors.CalcNumHBA(mol),
'NumRings': rdMolDescriptors.CalcNumRings(mol),
'NumAromaticRings': rdMolDescriptors.CalcNumAromaticRings(mol),
'TPSA': Descriptors.TPSA(mol),
'LogP': Descriptors.MolLogP(mol),
'FractionCSP3': rdMolDescriptors.CalcFractionCsp3(mol)
}
except:
return None
print("Obliczam deskryptory RDKit dla poprawnych struktur...")
descriptors = df_valid['SMILES'].apply(calculate_rdkit_descriptors)
df_descriptors = descriptors.apply(pd.Series)
# Dodanie do głównego DataFrame
df_valid = pd.concat([df_valid.reset_index(drop=True),
df_descriptors.reset_index(drop=True)], axis=1)
print(f"✅ Obliczono {len(df_descriptors.columns)} deskryptorów dla {len(df_valid)} struktur")
print("\n📊 Statystyki deskryptorów molekularnych:")
print(df_descriptors.describe())
# ============================================================
# 6. WIZUALIZACJA WŁAŚCIWOŚCI MOLEKULARNYCH
# ============================================================
print("\n" + "="*70)
print("📊 WIZUALIZACJA WŁAŚCIWOŚCI MOLEKULARNYCH")
print("="*70)
fig, axes = plt.subplots(2, 3, figsize=(18, 10))
# 1. Masa molowa
axes[0, 0].hist(df_descriptors['MolWt'], bins=50, color='#3498db',
edgecolor='black', alpha=0.7)
axes[0, 0].set_xlabel('Masa molowa (g/mol)', fontweight='bold')
axes[0, 0].set_ylabel('Liczba cząsteczek', fontweight='bold')
axes[0, 0].set_title('Rozkład masy molowej', fontweight='bold', fontsize=12)
axes[0, 0].grid(axis='y', alpha=0.3)
# 2. Liczba ciężkich atomów
axes[0, 1].hist(df_descriptors['NumHeavyAtoms'], bins=range(1, 11),
color='#e74c3c', edgecolor='black', alpha=0.7)
axes[0, 1].set_xlabel('Liczba ciężkich atomów', fontweight='bold')
axes[0, 1].set_ylabel('Liczba cząsteczek', fontweight='bold')
axes[0, 1].set_title('Rozkład liczby ciężkich atomów (max=9)', fontweight='bold', fontsize=12)
axes[0, 1].grid(axis='y', alpha=0.3)
# 3. LogP
axes[0, 2].hist(df_descriptors['LogP'], bins=50, color='#2ecc71',
edgecolor='black', alpha=0.7)
axes[0, 2].set_xlabel('LogP', fontweight='bold')
axes[0, 2].set_ylabel('Liczba cząsteczek', fontweight='bold')
axes[0, 2].set_title('Rozkład lipofilowości (LogP)', fontweight='bold', fontsize=12)
axes[0, 2].grid(axis='y', alpha=0.3)
# 4. TPSA
axes[1, 0].hist(df_descriptors['TPSA'], bins=50, color='#f39c12',
edgecolor='black', alpha=0.7)
axes[1, 0].set_xlabel('TPSA (Ų)', fontweight='bold')
axes[1, 0].set_ylabel('Liczba cząsteczek', fontweight='bold')
axes[1, 0].set_title('Rozkład polarnej powierzchni (TPSA)', fontweight='bold', fontsize=12)
axes[1, 0].grid(axis='y', alpha=0.3)
# 5. Liczba pierścieni
ring_counts = df_descriptors['NumRings'].value_counts().sort_index()
axes[1, 1].bar(ring_counts.index, ring_counts.values, color='#9b59b6',
edgecolor='black', alpha=0.7)
axes[1, 1].set_xlabel('Liczba pierścieni', fontweight='bold')
axes[1, 1].set_ylabel('Liczba cząsteczek', fontweight='bold')
axes[1, 1].set_title('Rozkład liczby pierścieni', fontweight='bold', fontsize=12)
axes[1, 1].grid(axis='y', alpha=0.3)
# 6. Frakcja sp3
axes[1, 2].hist(df_descriptors['FractionCSP3'], bins=50, color='#1abc9c',
edgecolor='black', alpha=0.7)
axes[1, 2].set_xlabel('Frakcja węgli sp³', fontweight='bold')
axes[1, 2].set_ylabel('Liczba cząsteczek', fontweight='bold')
axes[1, 2].set_title('Rozkład frakcji sp³ (saturacja)', fontweight='bold', fontsize=12)
axes[1, 2].grid(axis='y', alpha=0.3)
plt.suptitle('Właściwości molekularne datasetu QM9', fontsize=16, fontweight='bold', y=1.00)
plt.tight_layout()
plt.savefig('qm9_molecular_properties.png', dpi=300, bbox_inches='tight')
print("✅ Zapisano: qm9_molecular_properties.png")
plt.show()
print("\n" + "="*70)
print("✅ WALIDACJA I ANALIZA SMILES ZAKOŃCZONA!")
print("="*70)
======================================================================
WALIDACJA STRUKTUR SMILES - QM9 DATASET
======================================================================
📋 ETAP 1: PODSTAWOWA WALIDACJA
----------------------------------------------------------------------
Walidacja struktur SMILES...
✅ Poprawne SMILES: 1,000 (100.00%)
❌ Niepoprawne SMILES: 0 (0.00%)
✓ Kontynuuję analizę 1,000 poprawnych struktur...
======================================================================
🧪 ANALIZA SKŁADU CHEMICZNEGO
======================================================================
Analizuję skład atomowy cząsteczek...
📊 Typy atomów występujące w datasecie (5):
• C występuje w 1,000 cząsteczkach (100.00%)
• O występuje w 843 cząsteczkach (84.30%)
• N występuje w 622 cząsteczkach (62.20%)
• H występuje w 16 cząsteczkach ( 1.60%)
• F występuje w 14 cząsteczkach ( 1.40%)
📈 Statystyki liczby atomów:
C N O H F
count 1000.000000 1000.000000 1000.000000 1000.000000 1000.000000
mean 6.251000 1.093000 1.392000 0.017000 0.020000
std 1.282046 1.138265 0.875844 0.136857 0.188774
min 2.000000 0.000000 0.000000 0.000000 0.000000
25% 5.000000 0.000000 1.000000 0.000000 0.000000
50% 6.000000 1.000000 1.000000 0.000000 0.000000
75% 7.000000 2.000000 2.000000 0.000000 0.000000
max 9.000000 6.000000 4.000000 2.000000 3.000000
======================================================================
📊 WIZUALIZACJA SKŁADU ATOMOWEGO
======================================================================
✅ Zapisano: qm9_atomic_composition.png
====================================================================== ⚗️ OBLICZANIE DESKRYPTORÓW MOLEKULARNYCH ====================================================================== Obliczam deskryptory RDKit dla poprawnych struktur... ✅ Obliczono 0 deskryptorów dla 1000 struktur 📊 Statystyki deskryptorów molekularnych:
--------------------------------------------------------------------------- ValueError Traceback (most recent call last) Cell In[17], line 213 210 print(f"✅ Obliczono {len(df_descriptors.columns)} deskryptorów dla {len(df_valid)} struktur") 212 print("\n📊 Statystyki deskryptorów molekularnych:") --> 213 print(df_descriptors.describe()) 215 # ============================================================ 216 # 6. WIZUALIZACJA WŁAŚCIWOŚCI MOLEKULARNYCH 217 # ============================================================ 219 print("\n" + "="*70) File c:\Users\slast\miniconda3\envs\chem_eda\lib\site-packages\pandas\core\generic.py:12041, in NDFrame.describe(self, percentiles, include, exclude) 11799 @final 11800 def describe( 11801 self, (...) 11804 exclude=None, 11805 ) -> Self: 11806 """ 11807 Generate descriptive statistics. 11808 (...) 12039 max NaN 3.0 12040 """ > 12041 return describe_ndframe( 12042 obj=self, 12043 include=include, 12044 exclude=exclude, 12045 percentiles=percentiles, 12046 ).__finalize__(self, method="describe") File c:\Users\slast\miniconda3\envs\chem_eda\lib\site-packages\pandas\core\methods\describe.py:91, in describe_ndframe(obj, include, exclude, percentiles) 87 describer = SeriesDescriber( 88 obj=cast("Series", obj), 89 ) 90 else: ---> 91 describer = DataFrameDescriber( 92 obj=cast("DataFrame", obj), 93 include=include, 94 exclude=exclude, 95 ) 97 result = describer.describe(percentiles=percentiles) 98 return cast(NDFrameT, result) File c:\Users\slast\miniconda3\envs\chem_eda\lib\site-packages\pandas\core\methods\describe.py:162, in DataFrameDescriber.__init__(self, obj, include, exclude) 159 self.exclude = exclude 161 if obj.ndim == 2 and obj.columns.size == 0: --> 162 raise ValueError("Cannot describe a DataFrame without columns") 164 super().__init__(obj) ValueError: Cannot describe a DataFrame without columns
do naprawy potem¶
Rozkład train/valid/test - czy zbalansowany¶
# ============================================================
# KOMÓRKA 8: ANALIZA PODZIAŁU TRAIN/VALID/TEST
# ============================================================
print("="*70)
print("ANALIZA PODZIAŁU DATASETU - TRAIN/VALID/TEST")
print("="*70)
# ============================================================
# 1. PODSTAWOWE STATYSTYKI PODZIAŁU
# ============================================================
print("\n📊 ETAP 1: ROZKŁAD DANYCH WEDŁUG ZBIORÓW")
print("-" * 70)
# Liczebność poszczególnych zbiorów
split_counts = df['dataset_split'].value_counts().sort_index()
total_samples = len(df)
print("\n📈 Liczebność zbiorów:")
for split, count in split_counts.items():
pct = (count / total_samples) * 100
print(f" • {split:6s}: {count:7,} próbek ({pct:5.2f}%)")
print(f"\n RAZEM: {total_samples:7,} próbek (100.00%)")
# ============================================================
# 2. WIZUALIZACJA PODZIAŁU
# ============================================================
print("\n" + "="*70)
print("📊 WIZUALIZACJA PODZIAŁU")
print("="*70)
fig, axes = plt.subplots(1, 3, figsize=(18, 5))
# 1. Wykres kołowy
colors = ['#3498db', '#2ecc71', '#e74c3c']
wedges, texts, autotexts = axes[0].pie(split_counts.values,
labels=split_counts.index,
autopct='%1.1f%%',
colors=colors,
startangle=90,
textprops={'fontsize': 12, 'fontweight': 'bold'})
axes[0].set_title('Rozkład procentowy zbiorów', fontsize=14, fontweight='bold')
# 2. Wykres słupkowy
bars = axes[1].bar(split_counts.index, split_counts.values,
color=colors, edgecolor='black', linewidth=1.5, alpha=0.8)
axes[1].set_ylabel('Liczba próbek', fontsize=12, fontweight='bold')
axes[1].set_xlabel('Zbiór', fontsize=12, fontweight='bold')
axes[1].set_title('Liczebność zbiorów', fontsize=14, fontweight='bold')
axes[1].grid(axis='y', alpha=0.3)
# Dodanie wartości na słupkach
for bar, val in zip(bars, split_counts.values):
height = bar.get_height()
axes[1].text(bar.get_x() + bar.get_width()/2., height,
f'{val:,}', ha='center', va='bottom', fontweight='bold', fontsize=11)
# 3. Proporcje względem całości
proportions = (split_counts / total_samples * 100).values
axes[2].barh(split_counts.index, proportions, color=colors,
edgecolor='black', linewidth=1.5, alpha=0.8)
axes[2].set_xlabel('Procent całego datasetu (%)', fontsize=12, fontweight='bold')
axes[2].set_title('Proporcje zbiorów', fontsize=14, fontweight='bold')
axes[2].grid(axis='x', alpha=0.3)
# Dodanie wartości
for i, (split, val) in enumerate(zip(split_counts.index, proportions)):
axes[2].text(val + 0.5, i, f'{val:.1f}%',
va='center', fontweight='bold', fontsize=11)
plt.tight_layout()
plt.savefig('qm9_dataset_split.png', dpi=300, bbox_inches='tight')
print("✅ Zapisano: qm9_dataset_split.png")
plt.show()
# ============================================================
# 3. ANALIZA BALANSU WŁAŚCIWOŚCI MIĘDZY ZBIORAMI
# ============================================================
print("\n" + "="*70)
print("🔬 ANALIZA BALANSU WŁAŚCIWOŚCI KWANTOWYCH")
print("="*70)
# Wybór kluczowych właściwości do analizy
key_properties = ['HOMO', 'LUMO', 'gap', 'dipole_moment', 'polarizability',
'U0', 'H', 'G']
print("\n📊 Porównanie średnich wartości właściwości między zbiorami:")
print("-" * 70)
# Tabela porównawcza
comparison_data = []
for prop in key_properties:
if prop in df.columns:
train_mean = df[df['dataset_split'] == 'train'][prop].mean()
valid_mean = df[df['dataset_split'] == 'valid'][prop].mean()
test_mean = df[df['dataset_split'] == 'test'][prop].mean()
# Obliczenie różnic procentowych
max_diff = max(abs(train_mean - valid_mean),
abs(train_mean - test_mean),
abs(valid_mean - test_mean))
avg_mean = (train_mean + valid_mean + test_mean) / 3
diff_pct = (max_diff / abs(avg_mean) * 100) if avg_mean != 0 else 0
comparison_data.append({
'Właściwość': prop,
'Train': train_mean,
'Valid': valid_mean,
'Test': test_mean,
'Max Diff %': diff_pct
})
print(f"\n{prop:25s}")
print(f" Train: {train_mean:10.4f} | Valid: {valid_mean:10.4f} | Test: {test_mean:10.4f}")
print(f" Max różnica: {diff_pct:.2f}%")
df_comparison = pd.DataFrame(comparison_data)
# Ocena balansu
print("\n" + "="*70)
print("📋 OCENA BALANSU ZBIORÓW")
print("="*70)
max_imbalance = df_comparison['Max Diff %'].max()
print(f"\n✓ Maksymalna różnica między zbiorami: {max_imbalance:.2f}%")
if max_imbalance < 1:
balance_status = "🟢 DOSKONAŁY BALANS"
balance_desc = "Zbiory są bardzo dobrze zbalansowane"
elif max_imbalance < 3:
balance_status = "🟢 DOBRY BALANS"
balance_desc = "Zbiory są dobrze zbalansowane"
elif max_imbalance < 5:
balance_status = "🟡 AKCEPTOWALNY BALANS"
balance_desc = "Niewielkie różnice między zbiorami"
else:
balance_status = "🔴 SŁABY BALANS"
balance_desc = "Znaczące różnice między zbiorami"
print(f"\n{balance_status}")
print(f" {balance_desc}")
# ============================================================
# 4. WIZUALIZACJA ROZKŁADÓW WŁAŚCIWOŚCI
# ============================================================
print("\n" + "="*70)
print("📊 WIZUALIZACJA ROZKŁADÓW KLUCZOWYCH WŁAŚCIWOŚCI")
print("="*70)
# Wybór 6 najważniejszych właściwości
plot_properties = ['HOMO', 'LUMO', 'gap', 'dipole_moment', 'U0', 'H']
fig, axes = plt.subplots(2, 3, figsize=(18, 10))
axes = axes.flatten()
for i, prop in enumerate(plot_properties):
if prop in df.columns:
# Rozkłady dla każdego zbioru
for split, color in zip(['train', 'valid', 'test'], colors):
data = df[df['dataset_split'] == split][prop].dropna()
axes[i].hist(data, bins=50, alpha=0.5, label=split,
color=color, edgecolor='black', linewidth=0.5)
axes[i].set_xlabel(prop, fontsize=11, fontweight='bold')
axes[i].set_ylabel('Liczba próbek', fontsize=11, fontweight='bold')
axes[i].set_title(f'Rozkład {prop} w zbiorach', fontsize=12, fontweight='bold')
axes[i].legend()
axes[i].grid(axis='y', alpha=0.3)
plt.suptitle('Porównanie rozkładów właściwości między zbiorami',
fontsize=16, fontweight='bold', y=1.00)
plt.tight_layout()
plt.savefig('qm9_split_distributions.png', dpi=300, bbox_inches='tight')
print("✅ Zapisano: qm9_split_distributions.png")
plt.show()
# ============================================================
# 5. BOXPLOTY - PORÓWNANIE ZBIORÓW
# ============================================================
print("\n" + "="*70)
print("📦 BOXPLOTY - PORÓWNANIE MEDIANOWE")
print("="*70)
fig, axes = plt.subplots(2, 3, figsize=(18, 10))
axes = axes.flatten()
for i, prop in enumerate(plot_properties):
if prop in df.columns:
# Przygotowanie danych
data_to_plot = [df[df['dataset_split'] == split][prop].dropna()
for split in ['train', 'valid', 'test']]
bp = axes[i].boxplot(data_to_plot,
labels=['Train', 'Valid', 'Test'],
patch_artist=True,
widths=0.6)
# Kolorowanie boxów
for patch, color in zip(bp['boxes'], colors):
patch.set_facecolor(color)
patch.set_alpha(0.7)
axes[i].set_ylabel(prop, fontsize=11, fontweight='bold')
axes[i].set_title(f'Boxplot {prop}', fontsize=12, fontweight='bold')
axes[i].grid(axis='y', alpha=0.3)
plt.suptitle('Boxploty właściwości - porównanie zbiorów',
fontsize=16, fontweight='bold', y=1.00)
plt.tight_layout()
plt.savefig('qm9_split_boxplots.png', dpi=300, bbox_inches='tight')
print("✅ Zapisano: qm9_split_boxplots.png")
plt.show()
# ============================================================
# 6. PORÓWNANIE ODCHYLEŃ STANDARDOWYCH
# ============================================================
print("\n" + "="*70)
print("📊 PORÓWNANIE STATYSTYK OPISOWYCH")
print("="*70)
print("\nOdchylenia standardowe (im bardziej podobne, tym lepszy balans):")
print("-" * 70)
for prop in key_properties:
if prop in df.columns:
train_std = df[df['dataset_split'] == 'train'][prop].std()
valid_std = df[df['dataset_split'] == 'valid'][prop].std()
test_std = df[df['dataset_split'] == 'test'][prop].std()
# Obliczenie różnicy
max_std = max(train_std, valid_std, test_std)
min_std = min(train_std, valid_std, test_std)
diff_pct = ((max_std - min_std) / max_std * 100) if max_std > 0 else 0
print(f"\n{prop:20s}")
print(f" Train: {train_std:10.4f} | Valid: {valid_std:10.4f} | Test: {test_std:10.4f}")
print(f" Różnica: {diff_pct:.2f}%")
# ============================================================
# 7. PODSUMOWANIE
# ============================================================
print("\n" + "="*70)
print("📋 PODSUMOWANIE ANALIZY PODZIAŁU")
print("="*70)
similar_count = (df_comparison['Max Diff %'] < 3).sum()
total_props = len(df_comparison)
print(f"""
📊 ROZKŁAD LICZBOWY:
• Train: {split_counts['train']:,} ({split_counts['train']/total_samples*100:.1f}%)
• Valid: {split_counts['valid']:,} ({split_counts['valid']/total_samples*100:.1f}%)
• Test: {split_counts['test']:,} ({split_counts['test']/total_samples*100:.1f}%)
⚖️ BALANS WŁAŚCIWOŚCI:
• Status: {balance_status}
• Maksymalna różnica średnich: {max_imbalance:.2f}%
• Właściwości z dobrym balansem: {similar_count}/{total_props}
✅ WNIOSKI:
• {'Podział jest dobrze zbalansowany pod względem liczebności' if max(split_counts.values)/min(split_counts.values) < 1.2 else 'Podział ma lekką nierówność liczebności'}
• {'Rozkłady właściwości są spójne między zbiorami' if max_imbalance < 3 else 'Niektóre właściwości mogą różnić się między zbiorami'}
• {'Dataset jest gotowy do trenowania modeli ML' if max_imbalance < 5 else 'Rozważ dodatkowe sprawdzenie podziału'}
💡 REKOMENDACJE:
• {'Brak konieczności re-splittingu' if max_imbalance < 3 else 'Rozważ stratified split dla krytycznych właściwości'}
• {'Możesz bezpiecznie trenować modele na tym podziale' if similar_count/total_props > 0.7 else 'Monitoruj metryki na wszystkich zbiorach podczas treningu'}
""")
print("="*70)
print("✅ ANALIZA PODZIAŁU ZAKOŃCZONA!")
print("="*70)
====================================================================== ANALIZA PODZIAŁU DATASETU - TRAIN/VALID/TEST ====================================================================== 📊 ETAP 1: ROZKŁAD DANYCH WEDŁUG ZBIORÓW ---------------------------------------------------------------------- 📈 Liczebność zbiorów: • test : 91 próbek ( 9.10%) • train : 813 próbek (81.30%) • valid : 96 próbek ( 9.60%) RAZEM: 1,000 próbek (100.00%) ====================================================================== 📊 WIZUALIZACJA PODZIAŁU ====================================================================== ✅ Zapisano: qm9_dataset_split.png
====================================================================== 🔬 ANALIZA BALANSU WŁAŚCIWOŚCI KWANTOWYCH ====================================================================== 📊 Porównanie średnich wartości właściwości między zbiorami: ---------------------------------------------------------------------- HOMO Train: -0.0224 | Valid: -0.0890 | Test: -0.0519 Max różnica: 122.23% LUMO Train: 0.0727 | Valid: -0.1290 | Test: -0.1003 Max różnica: 386.37% gap Train: 0.0839 | Valid: -0.0843 | Test: -0.0742 Max różnica: 676.65% dipole_moment Train: -0.0299 | Valid: -0.0693 | Test: -0.0120 Max różnica: 154.55% polarizability Train: -0.0297 | Valid: -0.0061 | Test: -0.0796 Max różnica: 191.28% U0 Train: -0.0125 | Valid: -0.0063 | Test: -0.0707 Max różnica: 215.82% H Train: 0.0289 | Valid: 0.1364 | Test: -0.0878 Max różnica: 867.44% G Train: 0.0289 | Valid: 0.1364 | Test: -0.0878 Max różnica: 867.44% ====================================================================== 📋 OCENA BALANSU ZBIORÓW ====================================================================== ✓ Maksymalna różnica między zbiorami: 867.44% 🔴 SŁABY BALANS Znaczące różnice między zbiorami ====================================================================== 📊 WIZUALIZACJA ROZKŁADÓW KLUCZOWYCH WŁAŚCIWOŚCI ====================================================================== ✅ Zapisano: qm9_split_distributions.png
====================================================================== 📦 BOXPLOTY - PORÓWNANIE MEDIANOWE ======================================================================
C:\Users\slast\AppData\Local\Temp\ipykernel_6436\2569215989.py:200: MatplotlibDeprecationWarning: The 'labels' parameter of boxplot() has been renamed 'tick_labels' since Matplotlib 3.9; support for the old name will be dropped in 3.11. bp = axes[i].boxplot(data_to_plot, C:\Users\slast\AppData\Local\Temp\ipykernel_6436\2569215989.py:200: MatplotlibDeprecationWarning: The 'labels' parameter of boxplot() has been renamed 'tick_labels' since Matplotlib 3.9; support for the old name will be dropped in 3.11. bp = axes[i].boxplot(data_to_plot, C:\Users\slast\AppData\Local\Temp\ipykernel_6436\2569215989.py:200: MatplotlibDeprecationWarning: The 'labels' parameter of boxplot() has been renamed 'tick_labels' since Matplotlib 3.9; support for the old name will be dropped in 3.11. bp = axes[i].boxplot(data_to_plot, C:\Users\slast\AppData\Local\Temp\ipykernel_6436\2569215989.py:200: MatplotlibDeprecationWarning: The 'labels' parameter of boxplot() has been renamed 'tick_labels' since Matplotlib 3.9; support for the old name will be dropped in 3.11. bp = axes[i].boxplot(data_to_plot, C:\Users\slast\AppData\Local\Temp\ipykernel_6436\2569215989.py:200: MatplotlibDeprecationWarning: The 'labels' parameter of boxplot() has been renamed 'tick_labels' since Matplotlib 3.9; support for the old name will be dropped in 3.11. bp = axes[i].boxplot(data_to_plot, C:\Users\slast\AppData\Local\Temp\ipykernel_6436\2569215989.py:200: MatplotlibDeprecationWarning: The 'labels' parameter of boxplot() has been renamed 'tick_labels' since Matplotlib 3.9; support for the old name will be dropped in 3.11. bp = axes[i].boxplot(data_to_plot,
✅ Zapisano: qm9_split_boxplots.png
====================================================================== 📊 PORÓWNANIE STATYSTYK OPISOWYCH ====================================================================== Odchylenia standardowe (im bardziej podobne, tym lepszy balans): ---------------------------------------------------------------------- HOMO Train: 1.0495 | Valid: 0.9226 | Test: 0.9963 Różnica: 12.09% LUMO Train: 0.9760 | Valid: 1.0129 | Test: 1.1192 Różnica: 12.80% gap Train: 1.0070 | Valid: 0.9814 | Test: 1.0664 Różnica: 7.97% dipole_moment Train: 1.0284 | Valid: 0.9200 | Test: 0.9418 Różnica: 10.55% polarizability Train: 1.0579 | Valid: 1.1738 | Test: 1.0389 Różnica: 11.50% U0 Train: 1.0377 | Valid: 1.1266 | Test: 1.0621 Różnica: 7.89% H Train: 1.0520 | Valid: 1.1977 | Test: 1.0812 Różnica: 12.17% G Train: 1.0520 | Valid: 1.1977 | Test: 1.0812 Różnica: 12.17% ====================================================================== 📋 PODSUMOWANIE ANALIZY PODZIAŁU ====================================================================== 📊 ROZKŁAD LICZBOWY: • Train: 813 (81.3%) • Valid: 96 (9.6%) • Test: 91 (9.1%) ⚖️ BALANS WŁAŚCIWOŚCI: • Status: 🔴 SŁABY BALANS • Maksymalna różnica średnich: 867.44% • Właściwości z dobrym balansem: 0/8 ✅ WNIOSKI: • Podział ma lekką nierówność liczebności • Niektóre właściwości mogą różnić się między zbiorami • Rozważ dodatkowe sprawdzenie podziału 💡 REKOMENDACJE: • Rozważ stratified split dla krytycznych właściwości • Monitoruj metryki na wszystkich zbiorach podczas treningu ====================================================================== ✅ ANALIZA PODZIAŁU ZAKOŃCZONA! ======================================================================
3. ANALIZA SKŁADU CHEMICZNEGO 🧪¶
- Rozkład atomów (H, C, N, O, F) - ile każdego typu
- Masa molowa - rozkład, median, percentyle
- Liczba atomów ciężkich - histogram, statystyki
- Liczba atomów wodoru - analiza
- Wzory sumaryczne - najczęstsze składy
# ============================================================
# KOMÓRKA 9: ANALIZA SKŁADU CHEMICZNEGO
# ============================================================
print("="*70)
print("ANALIZA SKŁADU CHEMICZNEGO - QM9 DATASET")
print("="*70)
# ============================================================
# 1. ROZKŁAD ATOMÓW (H, C, N, O, F)
# ============================================================
print("\n🧪 ETAP 1: ROZKŁAD ATOMÓW W DATASECIE")
print("-" * 70)
# Funkcja do analizy składu
def get_atom_composition(smiles):
"""Zwraca słownik z liczbą atomów każdego typu"""
from rdkit import Chem
mol = Chem.MolFromSmiles(smiles)
if mol is None:
return None
atom_dict = {'H': 0, 'C': 0, 'N': 0, 'O': 0, 'F': 0}
for atom in mol.GetAtoms():
symbol = atom.GetSymbol()
if symbol in atom_dict:
atom_dict[symbol] += 1
# Dodaj wodory (mogą być jawne lub niejawne)
atom_dict['H'] = sum(1 for atom in mol.GetAtoms() if atom.GetSymbol() == 'H')
if atom_dict['H'] == 0: # Jeśli wodory są niejawne
atom_dict['H'] = sum(atom.GetTotalNumHs() for atom in mol.GetAtoms())
return atom_dict
# Obliczenie składu dla wszystkich cząsteczek
print("Analizuję skład atomowy cząsteczek...")
compositions = df['SMILES'].apply(get_atom_composition)
df_composition = compositions.apply(pd.Series)
# Dodanie do głównego DataFrame
df['H_count'] = df_composition['H']
df['C_count'] = df_composition['C']
df['N_count'] = df_composition['N']
df['O_count'] = df_composition['O']
df['F_count'] = df_composition['F']
print("\n📊 Statystyki rozkładu atomów:")
print("-" * 70)
atom_stats = df_composition.describe()
print(atom_stats)
# Całkowita liczba każdego atomu w datasecie
print("\n📈 Całkowita liczba atomów w datasecie:")
for atom in ['H', 'C', 'N', 'O', 'F']:
total = df[f'{atom}_count'].sum()
avg_per_mol = df[f'{atom}_count'].mean()
print(f" • {atom}: {total:,} atomów | Średnio {avg_per_mol:.2f} na cząsteczkę")
# ============================================================
# 2. WIZUALIZACJA ROZKŁADU ATOMÓW
# ============================================================
print("\n" + "="*70)
print("📊 WIZUALIZACJA ROZKŁADU ATOMÓW")
print("="*70)
fig, axes = plt.subplots(2, 3, figsize=(18, 10))
axes = axes.flatten()
atom_colors = {
'H': '#E8E8E8', # Jasny szary (wodór)
'C': '#909090', # Ciemny szary (węgiel)
'N': '#3050F8', # Niebieski (azot)
'O': '#FF0D0D', # Czerwony (tlen)
'F': '#90E050' # Zielony (fluor)
}
for i, atom in enumerate(['H', 'C', 'N', 'O', 'F']):
col = f'{atom}_count'
data = df[col].dropna()
axes[i].hist(data, bins=50, color=atom_colors[atom],
edgecolor='black', alpha=0.7, linewidth=0.5)
axes[i].set_xlabel(f'Liczba atomów {atom}', fontsize=11, fontweight='bold')
axes[i].set_ylabel('Liczba cząsteczek', fontsize=11, fontweight='bold')
axes[i].set_title(f'Rozkład atomów {atom}\nŚrednia: {data.mean():.2f} | Max: {data.max():.0f}',
fontsize=12, fontweight='bold')
axes[i].grid(axis='y', alpha=0.3)
# Dodanie linii mediany
median_val = data.median()
axes[i].axvline(median_val, color='red', linestyle='--', linewidth=2,
label=f'Mediana: {median_val:.1f}')
axes[i].legend()
# Wykres porównawczy - średnie liczby atomów
ax_comparison = axes[5]
atoms = ['H', 'C', 'N', 'O', 'F']
means = [df[f'{atom}_count'].mean() for atom in atoms]
colors_list = [atom_colors[atom] for atom in atoms]
bars = ax_comparison.bar(atoms, means, color=colors_list,
edgecolor='black', linewidth=1.5, alpha=0.8)
ax_comparison.set_ylabel('Średnia liczba atomów', fontsize=11, fontweight='bold')
ax_comparison.set_xlabel('Typ atomu', fontsize=11, fontweight='bold')
ax_comparison.set_title('Średni skład atomowy cząsteczek QM9',
fontsize=12, fontweight='bold')
ax_comparison.grid(axis='y', alpha=0.3)
for bar, val in zip(bars, means):
ax_comparison.text(bar.get_x() + bar.get_width()/2, val + 0.2,
f'{val:.2f}', ha='center', fontweight='bold', fontsize=11)
plt.suptitle('Rozkład atomów w datasecie QM9', fontsize=16, fontweight='bold', y=1.00)
plt.tight_layout()
plt.savefig('qm9_atom_distribution.png', dpi=300, bbox_inches='tight')
print("✅ Zapisano: qm9_atom_distribution.png")
plt.show()
# ============================================================
# 3. MASA MOLOWA - SZCZEGÓŁOWA ANALIZA
# ============================================================
print("\n" + "="*70)
print("⚖️ ANALIZA MASY MOLOWEJ")
print("="*70)
# Obliczenie masy molowej
from rdkit.Chem import Descriptors
def calculate_molecular_weight(smiles):
from rdkit import Chem
mol = Chem.MolFromSmiles(smiles)
if mol is None:
return None
return Descriptors.MolWt(mol)
df['MolWt'] = df['SMILES'].apply(calculate_molecular_weight)
# Statystyki
print("\n📊 Statystyki masy molowej:")
print(f" • Średnia: {df['MolWt'].mean():.2f} g/mol")
print(f" • Mediana: {df['MolWt'].median():.2f} g/mol")
print(f" • Odch. std: {df['MolWt'].std():.2f} g/mol")
print(f" • Min: {df['MolWt'].min():.2f} g/mol")
print(f" • Max: {df['MolWt'].max():.2f} g/mol")
# Percentyle
print("\n📈 Percentyle masy molowej:")
percentiles = [10, 25, 50, 75, 90, 95, 99]
for p in percentiles:
val = df['MolWt'].quantile(p/100)
print(f" • {p:2d}%: {val:.2f} g/mol")
# Wizualizacja
fig, axes = plt.subplots(1, 3, figsize=(18, 5))
# 1. Histogram
axes[0].hist(df['MolWt'], bins=60, color='steelblue',
edgecolor='black', alpha=0.7, linewidth=0.5)
axes[0].axvline(df['MolWt'].mean(), color='red', linestyle='--',
linewidth=2, label=f"Średnia: {df['MolWt'].mean():.1f}")
axes[0].axvline(df['MolWt'].median(), color='green', linestyle='--',
linewidth=2, label=f"Mediana: {df['MolWt'].median():.1f}")
axes[0].set_xlabel('Masa molowa (g/mol)', fontsize=12, fontweight='bold')
axes[0].set_ylabel('Liczba cząsteczek', fontsize=12, fontweight='bold')
axes[0].set_title('Rozkład masy molowej', fontsize=13, fontweight='bold')
axes[0].legend()
axes[0].grid(axis='y', alpha=0.3)
# 2. Boxplot
bp = axes[1].boxplot([df['MolWt'].dropna()], labels=['Masa molowa'],
patch_artist=True, widths=0.5)
bp['boxes'][0].set_facecolor('lightblue')
bp['boxes'][0].set_alpha(0.7)
axes[1].set_ylabel('Masa molowa (g/mol)', fontsize=12, fontweight='bold')
axes[1].set_title('Boxplot masy molowej', fontsize=13, fontweight='bold')
axes[1].grid(axis='y', alpha=0.3)
# 3. Dystrybuanta (CDF)
sorted_mw = np.sort(df['MolWt'].dropna())
cdf = np.arange(1, len(sorted_mw) + 1) / len(sorted_mw)
axes[2].plot(sorted_mw, cdf, linewidth=2, color='darkblue')
axes[2].set_xlabel('Masa molowa (g/mol)', fontsize=12, fontweight='bold')
axes[2].set_ylabel('Dystrybuanta (CDF)', fontsize=12, fontweight='bold')
axes[2].set_title('Dystrybuanta masy molowej', fontsize=13, fontweight='bold')
axes[2].grid(alpha=0.3)
# Dodanie linii percentyli
for p in [25, 50, 75]:
val = df['MolWt'].quantile(p/100)
axes[2].axvline(val, color='red', linestyle=':', alpha=0.5)
axes[2].text(val, 0.05, f'{p}%', rotation=90, fontsize=9)
plt.tight_layout()
plt.savefig('qm9_molecular_weight.png', dpi=300, bbox_inches='tight')
print("\n✅ Zapisano: qm9_molecular_weight.png")
plt.show()
# ============================================================
# 4. LICZBA ATOMÓW CIĘŻKICH
# ============================================================
print("\n" + "="*70)
print("🔬 ANALIZA ATOMÓW CIĘŻKICH (non-H)")
print("="*70)
# Obliczenie liczby atomów ciężkich
df['HeavyAtoms'] = df['C_count'] + df['N_count'] + df['O_count'] + df['F_count']
print("\n📊 Statystyki atomów ciężkich:")
print(f" • Średnia: {df['HeavyAtoms'].mean():.2f} atomów")
print(f" • Mediana: {df['HeavyAtoms'].median():.0f} atomów")
print(f" • Min: {df['HeavyAtoms'].min():.0f} atomów")
print(f" • Max: {df['HeavyAtoms'].max():.0f} atomów")
print(f" • Limit QM9: 9 atomów ciężkich")
# Rozkład
print("\n📈 Rozkład liczby atomów ciężkich:")
heavy_dist = df['HeavyAtoms'].value_counts().sort_index()
for count, freq in heavy_dist.items():
pct = (freq / len(df)) * 100
bar = '█' * int(pct / 2)
print(f" {count:2.0f} atomów: {freq:6,} cząsteczek ({pct:5.2f}%) {bar}")
# Wizualizacja
fig, axes = plt.subplots(1, 2, figsize=(14, 5))
# Histogram
axes[0].bar(heavy_dist.index, heavy_dist.values, color='coral',
edgecolor='black', alpha=0.8, linewidth=1)
axes[0].set_xlabel('Liczba atomów ciężkich', fontsize=12, fontweight='bold')
axes[0].set_ylabel('Liczba cząsteczek', fontsize=12, fontweight='bold')
axes[0].set_title('Rozkład liczby atomów ciężkich\n(QM9: max 9 atomów)',
fontsize=13, fontweight='bold')
axes[0].grid(axis='y', alpha=0.3)
axes[0].axvline(9, color='red', linestyle='--', linewidth=2, label='Limit QM9: 9')
axes[0].legend()
# Pie chart dla grup
bins = [0, 3, 6, 9, 12]
labels = ['1-3 atomy', '4-6 atomów', '7-9 atomów', '10+ atomów']
df['HeavyAtoms_bin'] = pd.cut(df['HeavyAtoms'], bins=bins, labels=labels[:len(bins)-1])
heavy_grouped = df['HeavyAtoms_bin'].value_counts()
axes[1].pie(heavy_grouped.values, labels=heavy_grouped.index, autopct='%1.1f%%',
colors=['#ff9999', '#66b3ff', '#99ff99'], startangle=90,
textprops={'fontsize': 11, 'fontweight': 'bold'})
axes[1].set_title('Grupy wielkości molekuł', fontsize=13, fontweight='bold')
plt.tight_layout()
plt.savefig('qm9_heavy_atoms.png', dpi=300, bbox_inches='tight')
print("\n✅ Zapisano: qm9_heavy_atoms.png")
plt.show()
# ============================================================
# 5. WZORY SUMARYCZNE - NAJCZĘSTSZE SKŁADY
# ============================================================
print("\n" + "="*70)
print("🧬 WZORY SUMARYCZNE - NAJCZĘSTSZE SKŁADY")
print("="*70)
# Tworzenie wzorów sumarycznych
def create_formula(row):
"""Tworzy wzór sumaryczny w formacie CxHyNzOaFb"""
parts = []
for atom in ['C', 'H', 'N', 'O', 'F']:
count = int(row[f'{atom}_count'])
if count > 0:
if count == 1:
parts.append(atom)
else:
parts.append(f"{atom}{count}")
return ''.join(parts) if parts else "Empty"
df['MolecularFormula'] = df.apply(create_formula, axis=1)
# Top 20 najczęstszych wzorów
print("\n🏆 TOP 20 najczęstszych wzorów sumarycznych:")
print("-" * 70)
top_formulas = df['MolecularFormula'].value_counts().head(20)
for i, (formula, count) in enumerate(top_formulas.items(), 1):
pct = (count / len(df)) * 100
print(f" {i:2d}. {formula:15s} {count:4,} cząsteczek ({pct:5.2f}%)")
# Wizualizacja top 15
fig, axes = plt.subplots(1, 2, figsize=(16, 6))
# Barplot
top15 = top_formulas.head(15)
axes[0].barh(range(len(top15)), top15.values, color='teal',
edgecolor='black', alpha=0.7)
axes[0].set_yticks(range(len(top15)))
axes[0].set_yticklabels(top15.index, fontsize=9)
axes[0].set_xlabel('Liczba cząsteczek', fontsize=12, fontweight='bold')
axes[0].set_title('Top 15 najczęstszych wzorów sumarycznych',
fontsize=13, fontweight='bold')
axes[0].grid(axis='x', alpha=0.3)
axes[0].invert_yaxis()
# Dodanie wartości
for i, val in enumerate(top15.values):
axes[0].text(val + 5, i, f'{val}', va='center', fontsize=9, fontweight='bold')
# Rozkład unikalności
unique_formulas = df['MolecularFormula'].nunique()
total_molecules = len(df)
categories = ['Unikalne wzory', 'Powtarzające się']
values = [unique_formulas, total_molecules - unique_formulas]
axes[1].pie(values, labels=categories, autopct='%1.1f%%',
colors=['#ff6b6b', '#4ecdc4'], startangle=90,
textprops={'fontsize': 12, 'fontweight': 'bold'})
axes[1].set_title(f'Unikalność wzorów sumarycznych\n({unique_formulas:,} unikalnych / {total_molecules:,} total)',
fontsize=13, fontweight='bold')
plt.tight_layout()
plt.savefig('qm9_molecular_formulas.png', dpi=300, bbox_inches='tight')
print("\n✅ Zapisano: qm9_molecular_formulas.png")
plt.show()
# ============================================================
# 6. PODSUMOWANIE SKŁADU CHEMICZNEGO
# ============================================================
print("\n" + "="*70)
print("📋 PODSUMOWANIE SKŁADU CHEMICZNEGO")
print("="*70)
print(f"""
🧪 SKŁAD ATOMOWY:
• Średnia H: {df['H_count'].mean():.2f} (wodór)
• Średnia C: {df['C_count'].mean():.2f} (węgiel)
• Średnia N: {df['N_count'].mean():.2f} (azot)
• Średnia O: {df['O_count'].mean():.2f} (tlen)
• Średnia F: {df['F_count'].mean():.2f} (fluor)
⚖️ MASA MOLOWA:
• Zakres: {df['MolWt'].min():.1f} - {df['MolWt'].max():.1f} g/mol
• Średnia: {df['MolWt'].mean():.1f} g/mol
• Mediana: {df['MolWt'].median():.1f} g/mol
🔬 ATOMY CIĘŻKIE:
• Zakres: {df['HeavyAtoms'].min():.0f} - {df['HeavyAtoms'].max():.0f} atomów
• Średnia: {df['HeavyAtoms'].mean():.2f} atomów
• Limit QM9: 9 atomów ciężkich (spec datasetu)
🧬 WZORY SUMARYCZNE:
• Unikalne wzory: {unique_formulas:,}
• Całkowita liczba molekuł: {total_molecules:,}
• Najczęstszy wzór: {top_formulas.index[0]} ({top_formulas.values[0]} molekuł)
✅ CHARAKTERYSTYKA DATASETU:
• QM9 zawiera małe cząsteczki organiczne (max 9 ciężkich atomów)
• Dominują związki z H, C, N, O, F (jak w specyfikacji)
• Duża różnorodność strukturalna ({unique_formulas:,} unikalnych wzorów)
• Masa molowa typowo poniżej 150 g/mol (small molecules)
""")
print("="*70)
print("✅ ANALIZA SKŁADU CHEMICZNEGO ZAKOŃCZONA!")
print("="*70)
======================================================================
ANALIZA SKŁADU CHEMICZNEGO - QM9 DATASET
======================================================================
🧪 ETAP 1: ROZKŁAD ATOMÓW W DATASECIE
----------------------------------------------------------------------
Analizuję skład atomowy cząsteczek...
📊 Statystyki rozkładu atomów:
----------------------------------------------------------------------
H C N O F
count 1000.000000 1000.000000 1000.000000 1000.000000 1000.000000
mean 9.322000 6.251000 1.093000 1.392000 0.020000
std 3.103215 1.282046 1.138265 0.875844 0.188774
min 1.000000 2.000000 0.000000 0.000000 0.000000
25% 7.000000 5.000000 0.000000 1.000000 0.000000
50% 9.000000 6.000000 1.000000 1.000000 0.000000
75% 12.000000 7.000000 2.000000 2.000000 0.000000
max 18.000000 9.000000 6.000000 4.000000 3.000000
📈 Całkowita liczba atomów w datasecie:
• H: 9,322 atomów | Średnio 9.32 na cząsteczkę
• C: 6,251 atomów | Średnio 6.25 na cząsteczkę
• N: 1,093 atomów | Średnio 1.09 na cząsteczkę
• O: 1,392 atomów | Średnio 1.39 na cząsteczkę
• F: 20 atomów | Średnio 0.02 na cząsteczkę
======================================================================
📊 WIZUALIZACJA ROZKŁADU ATOMÓW
======================================================================
✅ Zapisano: qm9_atom_distribution.png
====================================================================== ⚖️ ANALIZA MASY MOLOWEJ ====================================================================== 📊 Statystyki masy molowej: • Średnia: 122.53 g/mol • Mediana: 125.13 g/mol • Odch. std: 8.67 g/mol • Min: 59.07 g/mol • Max: 142.12 g/mol 📈 Percentyle masy molowej: • 10%: 112.11 g/mol • 25%: 121.10 g/mol • 50%: 125.13 g/mol • 75%: 128.13 g/mol • 90%: 130.08 g/mol • 95%: 131.13 g/mol • 99%: 132.20 g/mol
C:\Users\slast\AppData\Local\Temp\ipykernel_6436\3117465052.py:175: MatplotlibDeprecationWarning: The 'labels' parameter of boxplot() has been renamed 'tick_labels' since Matplotlib 3.9; support for the old name will be dropped in 3.11. bp = axes[1].boxplot([df['MolWt'].dropna()], labels=['Masa molowa'],
✅ Zapisano: qm9_molecular_weight.png
====================================================================== 🔬 ANALIZA ATOMÓW CIĘŻKICH (non-H) ====================================================================== 📊 Statystyki atomów ciężkich: • Średnia: 8.76 atomów • Mediana: 9 atomów • Min: 4 atomów • Max: 9 atomów • Limit QM9: 9 atomów ciężkich 📈 Rozkład liczby atomów ciężkich: 4 atomów: 1 cząsteczek ( 0.10%) 5 atomów: 2 cząsteczek ( 0.20%) 6 atomów: 8 cząsteczek ( 0.80%) 7 atomów: 30 cząsteczek ( 3.00%) █ 8 atomów: 147 cząsteczek (14.70%) ███████ 9 atomów: 812 cząsteczek (81.20%) ████████████████████████████████████████ ✅ Zapisano: qm9_heavy_atoms.png
====================================================================== 🧬 WZORY SUMARYCZNE - NAJCZĘSTSZE SKŁADY ====================================================================== 🏆 TOP 20 najczęstszych wzorów sumarycznych: ---------------------------------------------------------------------- 1. C7H10O2 52 cząsteczek ( 5.20%) 2. C6H9NO2 39 cząsteczek ( 3.90%) 3. C7H12O2 36 cząsteczek ( 3.60%) 4. C8H12O 35 cząsteczek ( 3.50%) 5. C7H9NO 25 cząsteczek ( 2.50%) 6. C6H8O3 23 cząsteczek ( 2.30%) 7. C8H10O 21 cząsteczek ( 2.10%) 8. C6H7NO2 19 cząsteczek ( 1.90%) 9. C7H8O2 19 cząsteczek ( 1.90%) 10. C7H14NO 18 cząsteczek ( 1.80%) 11. C6H10N2O 17 cząsteczek ( 1.70%) 12. C6H10O3 16 cząsteczek ( 1.60%) 13. C7H12NO 15 cząsteczek ( 1.50%) 14. C8H14O 14 cząsteczek ( 1.40%) 15. C6H8N2O 14 cząsteczek ( 1.40%) 16. C6H11NO2 13 cząsteczek ( 1.30%) 17. C6H12O3 13 cząsteczek ( 1.30%) 18. C8H16O 12 cząsteczek ( 1.20%) 19. C7H7NO 12 cząsteczek ( 1.20%) 20. C9H14 11 cząsteczek ( 1.10%) ✅ Zapisano: qm9_molecular_formulas.png
====================================================================== 📋 PODSUMOWANIE SKŁADU CHEMICZNEGO ====================================================================== 🧪 SKŁAD ATOMOWY: • Średnia H: 9.32 (wodór) • Średnia C: 6.25 (węgiel) • Średnia N: 1.09 (azot) • Średnia O: 1.39 (tlen) • Średnia F: 0.02 (fluor) ⚖️ MASA MOLOWA: • Zakres: 59.1 - 142.1 g/mol • Średnia: 122.5 g/mol • Mediana: 125.1 g/mol 🔬 ATOMY CIĘŻKIE: • Zakres: 4 - 9 atomów • Średnia: 8.76 atomów • Limit QM9: 9 atomów ciężkich (spec datasetu) 🧬 WZORY SUMARYCZNE: • Unikalne wzory: 243 • Całkowita liczba molekuł: 1,000 • Najczęstszy wzór: C7H10O2 (52 molekuł) ✅ CHARAKTERYSTYKA DATASETU: • QM9 zawiera małe cząsteczki organiczne (max 9 ciężkich atomów) • Dominują związki z H, C, N, O, F (jak w specyfikacji) • Duża różnorodność strukturalna (243 unikalnych wzorów) • Masa molowa typowo poniżej 150 g/mol (small molecules) ====================================================================== ✅ ANALIZA SKŁADU CHEMICZNEGO ZAKOŃCZONA! ======================================================================
4. ANALIZA STRUKTURALNA 🔬¶
- Liczba wiązań - single, double, triple, aromatic
- Liczba pierścieni - histogram, rozkład 0-3 pierścienie
- Pierścienie aromatyczne vs alifatyczne
- Rotowalne wiązania - flexibility molekuł
- Frakcja sp3 - saturacja vs aromatyczność
- Chiralność - czy są centra chiralne
# ============================================================
# KOMÓRKA 10: ANALIZA STRUKTURALNA
# ============================================================
print("="*70)
print("ANALIZA STRUKTURALNA - QM9 DATASET")
print("="*70)
# ============================================================
# 1. ANALIZA WIĄZAŃ (SINGLE, DOUBLE, TRIPLE, AROMATIC)
# ============================================================
print("\n🔗 ETAP 1: ANALIZA TYPÓW WIĄZAŃ")
print("-" * 70)
from rdkit import Chem
from rdkit.Chem import rdMolDescriptors, Descriptors
def analyze_bonds(smiles):
"""Analiza typów wiązań w cząsteczce"""
mol = Chem.MolFromSmiles(smiles)
if mol is None:
return None
bond_types = {
'single': 0,
'double': 0,
'triple': 0,
'aromatic': 0
}
for bond in mol.GetBonds():
bond_type = bond.GetBondType()
if bond_type == Chem.BondType.SINGLE:
if not bond.GetIsAromatic():
bond_types['single'] += 1
elif bond_type == Chem.BondType.DOUBLE:
if not bond.GetIsAromatic():
bond_types['double'] += 1
elif bond_type == Chem.BondType.TRIPLE:
bond_types['triple'] += 1
if bond.GetIsAromatic():
bond_types['aromatic'] += 1
return bond_types
print("Analizuję typy wiązań w cząsteczkach...")
bonds_data = df['SMILES'].apply(analyze_bonds)
df_bonds = bonds_data.apply(pd.Series)
# Dodanie do głównego DataFrame
df['single_bonds'] = df_bonds['single']
df['double_bonds'] = df_bonds['double']
df['triple_bonds'] = df_bonds['triple']
df['aromatic_bonds'] = df_bonds['aromatic']
print("\n📊 Statystyki typów wiązań:")
print(df_bonds.describe())
print("\n📈 Średnia liczba wiązań każdego typu:")
for bond_type in ['single', 'double', 'triple', 'aromatic']:
avg = df[f'{bond_type}_bonds'].mean()
total = df[f'{bond_type}_bonds'].sum()
molecules_with = (df[f'{bond_type}_bonds'] > 0).sum()
pct = (molecules_with / len(df)) * 100
print(f" • {bond_type:10s}: Średnio {avg:.2f} | Molekuły z tym wiązaniem: {molecules_with:,} ({pct:.1f}%)")
# Wizualizacja
fig, axes = plt.subplots(2, 2, figsize=(14, 10))
bond_colors = {
'single': '#3498db',
'double': '#e74c3c',
'triple': '#f39c12',
'aromatic': '#9b59b6'
}
for idx, (bond_type, ax) in enumerate(zip(['single', 'double', 'triple', 'aromatic'], axes.flatten())):
data = df[f'{bond_type}_bonds'].dropna()
ax.hist(data, bins=40, color=bond_colors[bond_type],
edgecolor='black', alpha=0.7, linewidth=0.5)
ax.set_xlabel(f'Liczba wiązań {bond_type}', fontsize=11, fontweight='bold')
ax.set_ylabel('Liczba cząsteczek', fontsize=11, fontweight='bold')
ax.set_title(f'Rozkład wiązań {bond_type}\nŚrednia: {data.mean():.2f}',
fontsize=12, fontweight='bold')
ax.grid(axis='y', alpha=0.3)
ax.axvline(data.mean(), color='red', linestyle='--', linewidth=2, label='Średnia')
ax.legend()
plt.suptitle('Rozkład typów wiązań w QM9', fontsize=16, fontweight='bold')
plt.tight_layout()
plt.savefig('qm9_bond_types.png', dpi=300, bbox_inches='tight')
print("\n✅ Zapisano: qm9_bond_types.png")
plt.show()
# ============================================================
# 2. LICZBA PIERŚCIENI
# ============================================================
print("\n" + "="*70)
print("⭕ ANALIZA PIERŚCIENI")
print("="*70)
def analyze_rings(smiles):
"""Szczegółowa analiza pierścieni"""
mol = Chem.MolFromSmiles(smiles)
if mol is None:
return None
return {
'total_rings': rdMolDescriptors.CalcNumRings(mol),
'aromatic_rings': rdMolDescriptors.CalcNumAromaticRings(mol),
'aliphatic_rings': rdMolDescriptors.CalcNumAliphaticRings(mol),
'saturated_rings': rdMolDescriptors.CalcNumSaturatedRings(mol),
'hetero_rings': rdMolDescriptors.CalcNumHeterocycles(mol),
'aromatic_hetero_rings': rdMolDescriptors.CalcNumAromaticHeterocycles(mol)
}
print("Analizuję struktury pierścieniowe...")
rings_data = df['SMILES'].apply(analyze_rings)
df_rings = rings_data.apply(pd.Series)
# Dodanie do głównego DataFrame
for col in df_rings.columns:
df[col] = df_rings[col]
print("\n📊 Statystyki pierścieni:")
print(df_rings.describe())
# Rozkład liczby pierścieni
print("\n📈 Rozkład liczby pierścieni:")
ring_dist = df['total_rings'].value_counts().sort_index()
for count, freq in ring_dist.items():
pct = (freq / len(df)) * 100
bar = '█' * int(pct / 3)
print(f" {count:1.0f} pierścieni: {freq:6,} cząsteczek ({pct:5.2f}%) {bar}")
# Molekuły z/bez pierścieni
molecules_with_rings = (df['total_rings'] > 0).sum()
molecules_without_rings = (df['total_rings'] == 0).sum()
print(f"\n✓ Z pierścieniami: {molecules_with_rings:,} ({molecules_with_rings/len(df)*100:.1f}%)")
print(f"✓ Bez pierścieni: {molecules_without_rings:,} ({molecules_without_rings/len(df)*100:.1f}%)")
# Wizualizacja
fig, axes = plt.subplots(2, 2, figsize=(14, 10))
# 1. Histogram liczby pierścieni
axes[0, 0].bar(ring_dist.index, ring_dist.values, color='steelblue',
edgecolor='black', alpha=0.8, linewidth=1)
axes[0, 0].set_xlabel('Liczba pierścieni', fontsize=12, fontweight='bold')
axes[0, 0].set_ylabel('Liczba cząsteczek', fontsize=12, fontweight='bold')
axes[0, 0].set_title('Rozkład liczby pierścieni', fontsize=13, fontweight='bold')
axes[0, 0].grid(axis='y', alpha=0.3)
# 2. Pie chart - z/bez pierścieni
axes[0, 1].pie([molecules_with_rings, molecules_without_rings],
labels=['Z pierścieniami', 'Bez pierścieni'],
autopct='%1.1f%%', colors=['#3498db', '#95a5a6'],
startangle=90, textprops={'fontsize': 11, 'fontweight': 'bold'})
axes[0, 1].set_title('Obecność pierścieni', fontsize=13, fontweight='bold')
# 3. Aromatyczne vs Alifatyczne
ring_types = ['aromatic_rings', 'aliphatic_rings', 'saturated_rings']
ring_labels = ['Aromatyczne', 'Alifatyczne', 'Nasycone']
ring_colors = ['#e74c3c', '#2ecc71', '#f39c12']
for ring_type, label, color in zip(ring_types, ring_labels, ring_colors):
data = df[ring_type].value_counts().sort_index()
axes[1, 0].plot(data.index, data.values, marker='o', linewidth=2,
label=label, color=color, markersize=8)
axes[1, 0].set_xlabel('Liczba pierścieni', fontsize=12, fontweight='bold')
axes[1, 0].set_ylabel('Liczba cząsteczek', fontsize=12, fontweight='bold')
axes[1, 0].set_title('Typy pierścieni', fontsize=13, fontweight='bold')
axes[1, 0].legend()
axes[1, 0].grid(alpha=0.3)
# 4. Heterocykle
hetero_dist = df['hetero_rings'].value_counts().sort_index()
axes[1, 1].bar(hetero_dist.index, hetero_dist.values, color='#9b59b6',
edgecolor='black', alpha=0.8, linewidth=1)
axes[1, 1].set_xlabel('Liczba heterocykli', fontsize=12, fontweight='bold')
axes[1, 1].set_ylabel('Liczba cząsteczek', fontsize=12, fontweight='bold')
axes[1, 1].set_title('Rozkład heterocykli', fontsize=13, fontweight='bold')
axes[1, 1].grid(axis='y', alpha=0.3)
plt.suptitle('Analiza struktur pierścieniowych', fontsize=16, fontweight='bold')
plt.tight_layout()
plt.savefig('qm9_rings_analysis.png', dpi=300, bbox_inches='tight')
print("\n✅ Zapisano: qm9_rings_analysis.png")
plt.show()
# ============================================================
# 3. ROTOWALNE WIĄZANIA (FLEXIBILITY)
# ============================================================
print("\n" + "="*70)
print("🔄 ANALIZA ROTOWALNYCH WIĄZAŃ (FLEXIBILITY)")
print("="*70)
def calculate_rotatable_bonds(smiles):
"""Oblicza liczbę rotowalnych wiązań"""
mol = Chem.MolFromSmiles(smiles)
if mol is None:
return None
return rdMolDescriptors.CalcNumRotatableBonds(mol)
df['rotatable_bonds'] = df['SMILES'].apply(calculate_rotatable_bonds)
print("\n📊 Statystyki rotowalnych wiązań:")
print(f" • Średnia: {df['rotatable_bonds'].mean():.2f} wiązań")
print(f" • Mediana: {df['rotatable_bonds'].median():.0f} wiązań")
print(f" • Max: {df['rotatable_bonds'].max():.0f} wiązań")
print(f" • Std: {df['rotatable_bonds'].std():.2f}")
# Rozkład
print("\n📈 Rozkład rotowalnych wiązań:")
rotatable_dist = df['rotatable_bonds'].value_counts().sort_index()
for count, freq in rotatable_dist.items():
pct = (freq / len(df)) * 100
if pct > 1: # Pokaż tylko istotne
bar = '█' * int(pct / 2)
print(f" {count:2.0f} wiązań: {freq:6,} cząsteczek ({pct:5.2f}%) {bar}")
# Klasyfikacja flexibility
rigid = (df['rotatable_bonds'] == 0).sum()
semi_flexible = ((df['rotatable_bonds'] > 0) & (df['rotatable_bonds'] <= 3)).sum()
flexible = (df['rotatable_bonds'] > 3).sum()
print(f"\n📋 Klasyfikacja elastyczności:")
print(f" • Sztywne (0 wiązań): {rigid:,} ({rigid/len(df)*100:.1f}%)")
print(f" • Semi-elastyczne (1-3): {semi_flexible:,} ({semi_flexible/len(df)*100:.1f}%)")
print(f" • Elastyczne (4+): {flexible:,} ({flexible/len(df)*100:.1f}%)")
# Wizualizacja
fig, axes = plt.subplots(1, 3, figsize=(18, 5))
# 1. Histogram
axes[0].hist(df['rotatable_bonds'], bins=range(0, int(df['rotatable_bonds'].max())+2),
color='teal', edgecolor='black', alpha=0.7)
axes[0].axvline(df['rotatable_bonds'].mean(), color='red', linestyle='--',
linewidth=2, label=f"Średnia: {df['rotatable_bonds'].mean():.2f}")
axes[0].set_xlabel('Liczba rotowalnych wiązań', fontsize=12, fontweight='bold')
axes[0].set_ylabel('Liczba cząsteczek', fontsize=12, fontweight='bold')
axes[0].set_title('Rozkład rotowalnych wiązań', fontsize=13, fontweight='bold')
axes[0].legend()
axes[0].grid(axis='y', alpha=0.3)
# 2. Pie chart - klasyfikacja
flexibility_labels = ['Sztywne\n(0)', 'Semi-elastyczne\n(1-3)', 'Elastyczne\n(4+)']
flexibility_values = [rigid, semi_flexible, flexible]
flexibility_colors = ['#e74c3c', '#f39c12', '#2ecc71']
axes[1].pie(flexibility_values, labels=flexibility_labels, autopct='%1.1f%%',
colors=flexibility_colors, startangle=90,
textprops={'fontsize': 11, 'fontweight': 'bold'})
axes[1].set_title('Klasyfikacja elastyczności molekuł', fontsize=13, fontweight='bold')
# 3. Zależność: rotatable bonds vs masa molowa
axes[2].scatter(df['MolWt'], df['rotatable_bonds'], alpha=0.3, s=10, color='darkblue')
axes[2].set_xlabel('Masa molowa (g/mol)', fontsize=12, fontweight='bold')
axes[2].set_ylabel('Liczba rotowalnych wiązań', fontsize=12, fontweight='bold')
axes[2].set_title('Elastyczność vs Masa molowa', fontsize=13, fontweight='bold')
axes[2].grid(alpha=0.3)
plt.tight_layout()
plt.savefig('qm9_rotatable_bonds.png', dpi=300, bbox_inches='tight')
print("\n✅ Zapisano: qm9_rotatable_bonds.png")
plt.show()
# ============================================================
# ============================================================
# 4. FRAKCJA SP3 (SATURACJA VS AROMATYCZNOŚĆ)
# ============================================================
print("\n" + "="*70)
print("🧬 ANALIZA FRAKCJI SP3 (SATURACJA)")
print("="*70)
def calculate_sp3_fraction(smiles):
"""Oblicza frakcję węgli sp3"""
mol = Chem.MolFromSmiles(smiles)
if mol is None:
return None
# Liczba węgli sp3
num_sp3 = 0
num_carbon = 0
for atom in mol.GetAtoms():
if atom.GetSymbol() == 'C':
num_carbon += 1
if atom.GetHybridization() == Chem.HybridizationType.SP3:
num_sp3 += 1
if num_carbon == 0:
return 0.0
return num_sp3 / num_carbon
df['frac_csp3'] = df['SMILES'].apply(calculate_sp3_fraction)
# ============================================================
# 5. CHIRALNOŚĆ
# ============================================================
print("\n" + "="*70)
print("🔀 ANALIZA CHIRALNOŚCI")
print("="*70)
def analyze_chirality(smiles):
"""Analiza centrów chiralnych"""
mol = Chem.MolFromSmiles(smiles)
if mol is None:
return None
# Znajdź centra chiralne
Chem.AssignStereochemistry(mol, cleanIt=True, force=True)
chiral_centers = Chem.FindMolChiralCenters(mol, includeUnassigned=True)
return {
'num_chiral_centers': len(chiral_centers),
'has_chirality': len(chiral_centers) > 0
}
print("Analizuję centra chiralne...")
chirality_data = df['SMILES'].apply(analyze_chirality)
df_chirality = chirality_data.apply(pd.Series)
df['num_chiral_centers'] = df_chirality['num_chiral_centers']
df['has_chirality'] = df_chirality['has_chirality']
# Statystyki
chiral_molecules = df['has_chirality'].sum()
achiral_molecules = (~df['has_chirality']).sum()
print(f"\n📊 Statystyki chiralności:")
print(f" • Cząsteczki chiralne: {chiral_molecules:,} ({chiral_molecules/len(df)*100:.1f}%)")
print(f" • Cząsteczki achiralne: {achiral_molecules:,} ({achiral_molecules/len(df)*100:.1f}%)")
print(f"\n📈 Rozkład liczby centrów chiralnych:")
chiral_dist = df['num_chiral_centers'].value_counts().sort_index()
for count, freq in chiral_dist.items():
pct = (freq / len(df)) * 100
if pct > 0.5:
bar = '█' * int(pct / 2)
print(f" {count:1.0f} centrów: {freq:6,} cząsteczek ({pct:5.2f}%) {bar}")
# Wizualizacja
fig, axes = plt.subplots(1, 3, figsize=(18, 5))
# 1. Pie chart - chiralne vs achiralne
axes[0].pie([chiral_molecules, achiral_molecules],
labels=['Chiralne', 'Achiralne'],
autopct='%1.1f%%', colors=['#e74c3c', '#95a5a6'],
startangle=90, textprops={'fontsize': 12, 'fontweight': 'bold'})
axes[0].set_title('Obecność chiralności', fontsize=13, fontweight='bold')
# 2. Histogram liczby centrów
chiral_only = df[df['has_chirality']]
if len(chiral_only) > 0:
axes[1].hist(chiral_only['num_chiral_centers'],
bins=range(0, int(chiral_only['num_chiral_centers'].max())+2),
color='coral', edgecolor='black', alpha=0.7)
axes[1].set_xlabel('Liczba centrów chiralnych', fontsize=12, fontweight='bold')
axes[1].set_ylabel('Liczba cząsteczek', fontsize=12, fontweight='bold')
axes[1].set_title('Rozkład centrów chiralnych\n(tylko molekuły chiralne)',
fontsize=13, fontweight='bold')
axes[1].grid(axis='y', alpha=0.3)
# 3. Chiralność vs złożoność
axes[2].scatter(df['HeavyAtoms'], df['num_chiral_centers'], alpha=0.3, s=10, color='darkred')
axes[2].set_xlabel('Liczba atomów ciężkich', fontsize=12, fontweight='bold')
axes[2].set_ylabel('Liczba centrów chiralnych', fontsize=12, fontweight='bold')
axes[2].set_title('Chiralność vs Złożoność molekuły', fontsize=13, fontweight='bold')
axes[2].grid(alpha=0.3)
plt.tight_layout()
plt.savefig('qm9_chirality.png', dpi=300, bbox_inches='tight')
print("\n✅ Zapisano: qm9_chirality.png")
plt.show()
# ============================================================
# 6. PODSUMOWANIE STRUKTURALNE
# ============================================================
print("\n" + "="*70)
print("📋 PODSUMOWANIE ANALIZY STRUKTURALNEJ")
print("="*70)
print(f"""
🔗 WIĄZANIA:
• Single: Średnio {df['single_bonds'].mean():.2f} wiązań
• Double: Średnio {df['double_bonds'].mean():.2f} wiązań
• Triple: Średnio {df['triple_bonds'].mean():.2f} wiązań
• Aromatic: Średnio {df['aromatic_bonds'].mean():.2f} wiązań
⭕ PIERŚCIENIE:
• Średnia liczba pierścieni: {df['total_rings'].mean():.2f}
• Molekuły z pierścieniami: {molecules_with_rings:,} ({molecules_with_rings/len(df)*100:.1f}%)
• Średnia pierścieni aromatycznych: {df['aromatic_rings'].mean():.2f}
• Średnia heterocykli: {df['hetero_rings'].mean():.2f}
🔄 FLEXIBILITY:
• Średnia rotowalnych wiązań: {df['rotatable_bonds'].mean():.2f}
• Molekuły sztywne (0): {rigid:,} ({rigid/len(df)*100:.1f}%)
• Molekuły elastyczne (4+): {flexible:,} ({flexible/len(df)*100:.1f}%)
🧬 NASYCENIE:
• Średnia frakcja sp³: {df['frac_csp3'].mean():.3f}
• Min frakcja sp³: {df['frac_csp3'].min():.3f}
• Max frakcja sp³: {df['frac_csp3'].max():.3f}
🔀 CHIRALNOŚĆ:
• Molekuły chiralne: {chiral_molecules:,} ({chiral_molecules/len(df)*100:.1f}%)
• Średnia centrów chiralnych: {df['num_chiral_centers'].mean():.2f}
✅ CHARAKTERYSTYKA:
• QM9 zawiera głównie małe, względnie sztywne molekuły
• Dominują struktury z 0-2 pierścieniami
• Mix molekuł aromatycznych i alifatycznych
• Niska chiralność (typowe dla małych molekuł)
""")
print("="*70)
print("✅ ANALIZA STRUKTURALNA ZAKOŃCZONA!")
print("="*70)
======================================================================
ANALIZA STRUKTURALNA - QM9 DATASET
======================================================================
🔗 ETAP 1: ANALIZA TYPÓW WIĄZAŃ
----------------------------------------------------------------------
Analizuję typy wiązań w cząsteczkach...
📊 Statystyki typów wiązań:
single double triple aromatic
count 1000.000000 1000.000000 1000.000000 1000.000000
mean 7.502000 0.613000 0.268000 0.974000
std 2.624287 0.749528 0.548157 2.129756
min 0.000000 0.000000 0.000000 0.000000
25% 6.000000 0.000000 0.000000 0.000000
50% 8.000000 0.000000 0.000000 0.000000
75% 9.000000 1.000000 0.000000 0.000000
max 13.000000 3.000000 4.000000 10.000000
📈 Średnia liczba wiązań każdego typu:
• single : Średnio 7.50 | Molekuły z tym wiązaniem: 997 (99.7%)
• double : Średnio 0.61 | Molekuły z tym wiązaniem: 472 (47.2%)
• triple : Średnio 0.27 | Molekuły z tym wiązaniem: 222 (22.2%)
• aromatic : Średnio 0.97 | Molekuły z tym wiązaniem: 179 (17.9%)
✅ Zapisano: qm9_bond_types.png
======================================================================
⭕ ANALIZA PIERŚCIENI
======================================================================
Analizuję struktury pierścieniowe...
📊 Statystyki pierścieni:
total_rings aromatic_rings aliphatic_rings saturated_rings \
count 1000.000000 1000.000000 1000.000000 1000.000000
mean 1.696000 0.189000 1.507000 1.333000
std 1.153653 0.416476 1.282028 1.302231
min 0.000000 0.000000 0.000000 0.000000
25% 1.000000 0.000000 1.000000 0.000000
50% 1.000000 0.000000 1.000000 1.000000
75% 2.000000 0.000000 2.000000 2.000000
max 7.000000 2.000000 7.000000 7.000000
hetero_rings aromatic_hetero_rings
count 1000.000000 1000.000000
mean 1.018000 0.188000
std 0.893574 0.415727
min 0.000000 0.000000
25% 0.000000 0.000000
50% 1.000000 0.000000
75% 1.000000 0.000000
max 6.000000 2.000000
📈 Rozkład liczby pierścieni:
0 pierścieni: 98 cząsteczek ( 9.80%) ███
1 pierścieni: 420 cząsteczek (42.00%) ██████████████
2 pierścieni: 271 cząsteczek (27.10%) █████████
3 pierścieni: 146 cząsteczek (14.60%) ████
4 pierścieni: 38 cząsteczek ( 3.80%) █
5 pierścieni: 21 cząsteczek ( 2.10%)
6 pierścieni: 3 cząsteczek ( 0.30%)
7 pierścieni: 3 cząsteczek ( 0.30%)
✓ Z pierścieniami: 902 (90.2%)
✓ Bez pierścieni: 98 (9.8%)
✅ Zapisano: qm9_rings_analysis.png
====================================================================== 🔄 ANALIZA ROTOWALNYCH WIĄZAŃ (FLEXIBILITY) ====================================================================== 📊 Statystyki rotowalnych wiązań: • Średnia: 0.89 wiązań • Mediana: 1 wiązań • Max: 5 wiązań • Std: 1.07 📈 Rozkład rotowalnych wiązań: 0 wiązań: 475 cząsteczek (47.50%) ███████████████████████ 1 wiązań: 281 cząsteczek (28.10%) ██████████████ 2 wiązań: 159 cząsteczek (15.90%) ███████ 3 wiązań: 55 cząsteczek ( 5.50%) ██ 4 wiązań: 24 cząsteczek ( 2.40%) █ 📋 Klasyfikacja elastyczności: • Sztywne (0 wiązań): 475 (47.5%) • Semi-elastyczne (1-3): 495 (49.5%) • Elastyczne (4+): 30 (3.0%) ✅ Zapisano: qm9_rotatable_bonds.png
====================================================================== 🧬 ANALIZA FRAKCJI SP3 (SATURACJA) ====================================================================== ====================================================================== 🔀 ANALIZA CHIRALNOŚCI ====================================================================== Analizuję centra chiralne... 📊 Statystyki chiralności: • Cząsteczki chiralne: 696 (69.6%) • Cząsteczki achiralne: 304 (30.4%) 📈 Rozkład liczby centrów chiralnych: 0 centrów: 304 cząsteczek (30.40%) ███████████████ 1 centrów: 184 cząsteczek (18.40%) █████████ 2 centrów: 203 cząsteczek (20.30%) ██████████ 3 centrów: 133 cząsteczek (13.30%) ██████ 4 centrów: 111 cząsteczek (11.10%) █████ 5 centrów: 44 cząsteczek ( 4.40%) ██ 6 centrów: 17 cząsteczek ( 1.70%) ✅ Zapisano: qm9_chirality.png
====================================================================== 📋 PODSUMOWANIE ANALIZY STRUKTURALNEJ ====================================================================== 🔗 WIĄZANIA: • Single: Średnio 7.50 wiązań • Double: Średnio 0.61 wiązań • Triple: Średnio 0.27 wiązań • Aromatic: Średnio 0.97 wiązań ⭕ PIERŚCIENIE: • Średnia liczba pierścieni: 1.70 • Molekuły z pierścieniami: 902 (90.2%) • Średnia pierścieni aromatycznych: 0.19 • Średnia heterocykli: 1.02 🔄 FLEXIBILITY: • Średnia rotowalnych wiązań: 0.89 • Molekuły sztywne (0): 475 (47.5%) • Molekuły elastyczne (4+): 30 (3.0%) 🧬 NASYCENIE: • Średnia frakcja sp³: 0.695 • Min frakcja sp³: 0.000 • Max frakcja sp³: 1.000 🔀 CHIRALNOŚĆ: • Molekuły chiralne: 696 (69.6%) • Średnia centrów chiralnych: 1.78 ✅ CHARAKTERYSTYKA: • QM9 zawiera głównie małe, względnie sztywne molekuły • Dominują struktury z 0-2 pierścieniami • Mix molekuł aromatycznych i alifatycznych • Niska chiralność (typowe dla małych molekuł) ====================================================================== ✅ ANALIZA STRUKTURALNA ZAKOŃCZONA! ======================================================================
5. WŁAŚCIWOŚCI FIZYKOCHEMICZNE (RDKit) ⚗️¶
- LogP - lipofilowość, rozkład, outliers
- TPSA - polar surface area, drug-likeness
- HBD/HBA - donory i akceptory wodoru
- Molar Refractivity - objętość molekularna
- Reguła Lipińskiego - ile struktur spełnia Ro5
- Korelacje między deskryptorami - heatmap
# ============================================================
# KOMÓRKA 11: WŁAŚCIWOŚCI FIZYKOCHEMICZNE (RDKit)
# ============================================================
print("="*70)
print("WŁAŚCIWOŚCI FIZYKOCHEMICZNE - QM9 DATASET")
print("="*70)
# ============================================================
# 1. OBLICZANIE DESKRYPTORÓW RDKIT
# ============================================================
print("\n⚗️ ETAP 1: OBLICZANIE DESKRYPTORÓW FIZYKOCHEMICZNYCH")
print("-" * 70)
from rdkit import Chem
from rdkit.Chem import Descriptors, Crippen, Lipinski
def calculate_physicochemical_properties(smiles):
"""Oblicza wszystkie kluczowe właściwości fizykochemiczne"""
mol = Chem.MolFromSmiles(smiles)
if mol is None:
return None
try:
return {
'LogP': Crippen.MolLogP(mol),
'MolMR': Crippen.MolMR(mol), # Molar Refractivity
'TPSA': Descriptors.TPSA(mol),
'HBD': Lipinski.NumHDonors(mol), # H-Bond Donors
'HBA': Lipinski.NumHAcceptors(mol), # H-Bond Acceptors
'LabuteASA': Descriptors.LabuteASA(mol), # Surface Area
'BalabanJ': Descriptors.BalabanJ(mol), # Topological index
}
except:
return None
print("Obliczam właściwości fizykochemiczne...")
props_data = df['SMILES'].apply(calculate_physicochemical_properties)
df_props = props_data.apply(pd.Series)
# Dodanie do głównego DataFrame
for col in df_props.columns:
df[col] = df_props[col]
print(f"\n✅ Obliczono {len(df_props.columns)} właściwości dla {len(df)} molekuł")
print("\n📊 Statystyki właściwości fizykochemicznych:")
print(df_props.describe())
# ============================================================
# 2. ANALIZA LogP (LIPOFILOWOŚĆ)
# ============================================================
print("\n" + "="*70)
print("💧 ANALIZA LogP (LIPOFILOWOŚĆ)")
print("="*70)
print("\n📊 Statystyki LogP:")
print(f" • Średnia: {df['LogP'].mean():.3f}")
print(f" • Mediana: {df['LogP'].median():.3f}")
print(f" • Std: {df['LogP'].std():.3f}")
print(f" • Min: {df['LogP'].min():.3f}")
print(f" • Max: {df['LogP'].max():.3f}")
# Percentyle
print("\n📈 Percentyle LogP:")
for p in [10, 25, 50, 75, 90]:
val = df['LogP'].quantile(p/100)
print(f" • {p:2d}%: {val:.3f}")
# Klasyfikacja lipofilowości
very_hydrophilic = (df['LogP'] < -1).sum()
hydrophilic = ((df['LogP'] >= -1) & (df['LogP'] < 1)).sum()
moderate = ((df['LogP'] >= 1) & (df['LogP'] < 3)).sum()
lipophilic = ((df['LogP'] >= 3) & (df['LogP'] < 5)).sum()
very_lipophilic = (df['LogP'] >= 5).sum()
print(f"\n📋 Klasyfikacja lipofilowości:")
print(f" • Bardzo hydrofilowe (LogP < -1): {very_hydrophilic:,} ({very_hydrophilic/len(df)*100:.1f}%)")
print(f" • Hydrofilowe (-1 ≤ LogP < 1): {hydrophilic:,} ({hydrophilic/len(df)*100:.1f}%)")
print(f" • Umiarkowane (1 ≤ LogP < 3): {moderate:,} ({moderate/len(df)*100:.1f}%)")
print(f" • Lipofilowe (3 ≤ LogP < 5): {lipophilic:,} ({lipophilic/len(df)*100:.1f}%)")
print(f" • Bardzo lipofilowe (LogP ≥ 5): {very_lipophilic:,} ({very_lipophilic/len(df)*100:.1f}%)")
# Outliers (IQR method)
Q1 = df['LogP'].quantile(0.25)
Q3 = df['LogP'].quantile(0.75)
IQR = Q3 - Q1
lower_bound = Q1 - 1.5 * IQR
upper_bound = Q3 + 1.5 * IQR
outliers_logp = df[(df['LogP'] < lower_bound) | (df['LogP'] > upper_bound)]
print(f"\n⚠️ Outliers LogP (metoda IQR):")
print(f" • Liczba outliers: {len(outliers_logp)} ({len(outliers_logp)/len(df)*100:.2f}%)")
print(f" • Zakres normalny: [{lower_bound:.2f}, {upper_bound:.2f}]")
# Wizualizacja
fig, axes = plt.subplots(2, 2, figsize=(14, 10))
# 1. Histogram
axes[0, 0].hist(df['LogP'], bins=60, color='dodgerblue', edgecolor='black', alpha=0.7)
axes[0, 0].axvline(df['LogP'].mean(), color='red', linestyle='--', linewidth=2,
label=f"Średnia: {df['LogP'].mean():.2f}")
axes[0, 0].axvline(df['LogP'].median(), color='green', linestyle='--', linewidth=2,
label=f"Mediana: {df['LogP'].median():.2f}")
axes[0, 0].set_xlabel('LogP', fontsize=12, fontweight='bold')
axes[0, 0].set_ylabel('Liczba molekuł', fontsize=12, fontweight='bold')
axes[0, 0].set_title('Rozkład LogP', fontsize=13, fontweight='bold')
axes[0, 0].legend()
axes[0, 0].grid(axis='y', alpha=0.3)
# 2. Boxplot z outliers
bp = axes[0, 1].boxplot([df['LogP'].dropna()], labels=['LogP'],
patch_artist=True, widths=0.5, showfliers=True)
bp['boxes'][0].set_facecolor('lightblue')
axes[0, 1].set_ylabel('LogP', fontsize=12, fontweight='bold')
axes[0, 1].set_title('Boxplot LogP (z outliers)', fontsize=13, fontweight='bold')
axes[0, 1].grid(axis='y', alpha=0.3)
# 3. Klasyfikacja - barplot
categories = ['Bardzo\nhydrofilowe', 'Hydrofilowe', 'Umiarkowane',
'Lipofilowe', 'Bardzo\nlipofilowe']
values = [very_hydrophilic, hydrophilic, moderate, lipophilic, very_lipophilic]
colors_logp = ['#3498db', '#5dade2', '#aed6f1', '#f8b739', '#e67e22']
bars = axes[1, 0].bar(range(len(categories)), values, color=colors_logp,
edgecolor='black', alpha=0.8, linewidth=1)
axes[1, 0].set_xticks(range(len(categories)))
axes[1, 0].set_xticklabels(categories, fontsize=9)
axes[1, 0].set_ylabel('Liczba molekuł', fontsize=12, fontweight='bold')
axes[1, 0].set_title('Klasyfikacja lipofilowości', fontsize=13, fontweight='bold')
axes[1, 0].grid(axis='y', alpha=0.3)
for bar, val in zip(bars, values):
axes[1, 0].text(bar.get_x() + bar.get_width()/2, val + 10,
f'{val}', ha='center', fontsize=9, fontweight='bold')
# 4. LogP vs Masa molowa
axes[1, 1].scatter(df['MolWt'], df['LogP'], alpha=0.3, s=10, color='darkblue')
axes[1, 1].set_xlabel('Masa molowa (g/mol)', fontsize=12, fontweight='bold')
axes[1, 1].set_ylabel('LogP', fontsize=12, fontweight='bold')
axes[1, 1].set_title('LogP vs Masa molowa', fontsize=13, fontweight='bold')
axes[1, 1].grid(alpha=0.3)
plt.suptitle('Analiza lipofilowości (LogP)', fontsize=16, fontweight='bold')
plt.tight_layout()
plt.savefig('qm9_logp_analysis.png', dpi=300, bbox_inches='tight')
print("\n✅ Zapisano: qm9_logp_analysis.png")
plt.show()
# ============================================================
# 3. ANALIZA TPSA (POLAR SURFACE AREA)
# ============================================================
print("\n" + "="*70)
print("🧊 ANALIZA TPSA (TOPOLOGICAL POLAR SURFACE AREA)")
print("="*70)
print("\n📊 Statystyki TPSA:")
print(f" • Średnia: {df['TPSA'].mean():.2f} Ų")
print(f" • Mediana: {df['TPSA'].median():.2f} Ų")
print(f" • Std: {df['TPSA'].std():.2f} Ų")
print(f" • Max: {df['TPSA'].max():.2f} Ų")
# Drug-likeness według TPSA
# TPSA < 140 Ų - dobra przepuszczalność BBB (blood-brain barrier)
# TPSA < 90 Ų - optymalna dla leków CNS
bbb_permeable = (df['TPSA'] < 140).sum()
cns_optimal = (df['TPSA'] < 90).sum()
print(f"\n💊 Drug-likeness (TPSA):")
print(f" • BBB permeable (< 140 Ų): {bbb_permeable:,} ({bbb_permeable/len(df)*100:.1f}%)")
print(f" • CNS optimal (< 90 Ų): {cns_optimal:,} ({cns_optimal/len(df)*100:.1f}%)")
# Wizualizacja
fig, axes = plt.subplots(1, 3, figsize=(18, 5))
# 1. Histogram
axes[0].hist(df['TPSA'], bins=50, color='coral', edgecolor='black', alpha=0.7)
axes[0].axvline(90, color='green', linestyle='--', linewidth=2, label='CNS optimal (90)')
axes[0].axvline(140, color='orange', linestyle='--', linewidth=2, label='BBB limit (140)')
axes[0].set_xlabel('TPSA (Ų)', fontsize=12, fontweight='bold')
axes[0].set_ylabel('Liczba molekuł', fontsize=12, fontweight='bold')
axes[0].set_title('Rozkład TPSA', fontsize=13, fontweight='bold')
axes[0].legend()
axes[0].grid(axis='y', alpha=0.3)
# 2. TPSA vs LogP
axes[1].scatter(df['LogP'], df['TPSA'], alpha=0.3, s=10, color='purple')
axes[1].set_xlabel('LogP', fontsize=12, fontweight='bold')
axes[1].set_ylabel('TPSA (Ų)', fontsize=12, fontweight='bold')
axes[1].set_title('TPSA vs LogP', fontsize=13, fontweight='bold')
axes[1].grid(alpha=0.3)
# 3. Pie chart - drug-likeness
categories_tpsa = ['CNS optimal\n(<90)', 'BBB permeable\n(90-140)', 'High TPSA\n(>140)']
values_tpsa = [cns_optimal, bbb_permeable - cns_optimal, len(df) - bbb_permeable]
colors_tpsa = ['#2ecc71', '#f39c12', '#e74c3c']
axes[2].pie(values_tpsa, labels=categories_tpsa, autopct='%1.1f%%',
colors=colors_tpsa, startangle=90,
textprops={'fontsize': 11, 'fontweight': 'bold'})
axes[2].set_title('TPSA Drug-likeness', fontsize=13, fontweight='bold')
plt.tight_layout()
plt.savefig('qm9_tpsa_analysis.png', dpi=300, bbox_inches='tight')
print("\n✅ Zapisano: qm9_tpsa_analysis.png")
plt.show()
# ============================================================
# 4. WIĄZANIA WODOROWE (HBD/HBA)
# ============================================================
print("\n" + "="*70)
print("🔗 ANALIZA WIĄZAŃ WODOROWYCH (HBD/HBA)")
print("="*70)
print("\n📊 Statystyki:")
print(f" • Średnia HBD (donory): {df['HBD'].mean():.2f}")
print(f" • Średnia HBA (akceptory): {df['HBA'].mean():.2f}")
print(f" • Max HBD: {df['HBD'].max():.0f}")
print(f" • Max HBA: {df['HBA'].max():.0f}")
# Rozkład
print("\n📈 Rozkład HBD:")
hbd_dist = df['HBD'].value_counts().sort_index()
for count, freq in hbd_dist.items():
pct = (freq / len(df)) * 100
if pct > 1:
bar = '█' * int(pct / 3)
print(f" {count:.0f} donorów: {freq:5,} ({pct:5.2f}%) {bar}")
print("\n📈 Rozkład HBA:")
hba_dist = df['HBA'].value_counts().sort_index()
for count, freq in hba_dist.items():
pct = (freq / len(df)) * 100
if pct > 1:
bar = '█' * int(pct / 3)
print(f" {count:.0f} akceptorów: {freq:5,} ({pct:5.2f}%) {bar}")
# Wizualizacja
fig, axes = plt.subplots(1, 3, figsize=(18, 5))
# 1. HBD histogram
axes[0].bar(hbd_dist.index, hbd_dist.values, color='#3498db',
edgecolor='black', alpha=0.8, linewidth=1)
axes[0].set_xlabel('Liczba donorów H (HBD)', fontsize=12, fontweight='bold')
axes[0].set_ylabel('Liczba molekuł', fontsize=12, fontweight='bold')
axes[0].set_title('Rozkład donorów wodoru', fontsize=13, fontweight='bold')
axes[0].grid(axis='y', alpha=0.3)
# 2. HBA histogram
axes[1].bar(hba_dist.index, hba_dist.values, color='#e74c3c',
edgecolor='black', alpha=0.8, linewidth=1)
axes[1].set_xlabel('Liczba akceptorów H (HBA)', fontsize=12, fontweight='bold')
axes[1].set_ylabel('Liczba molekuł', fontsize=12, fontweight='bold')
axes[1].set_title('Rozkład akceptorów wodoru', fontsize=13, fontweight='bold')
axes[1].grid(axis='y', alpha=0.3)
# 3. HBD vs HBA scatter
axes[2].scatter(df['HBD'], df['HBA'], alpha=0.4, s=15, color='green', edgecolors='black', linewidth=0.3)
axes[2].set_xlabel('HBD (donory)', fontsize=12, fontweight='bold')
axes[2].set_ylabel('HBA (akceptory)', fontsize=12, fontweight='bold')
axes[2].set_title('HBD vs HBA', fontsize=13, fontweight='bold')
axes[2].grid(alpha=0.3)
plt.tight_layout()
plt.savefig('qm9_hbond_analysis.png', dpi=300, bbox_inches='tight')
print("\n✅ Zapisano: qm9_hbond_analysis.png")
plt.show()
# ============================================================
# 5. REGUŁA LIPIŃSKIEGO (RULE OF 5)
# ============================================================
print("\n" + "="*70)
print("💊 REGUŁA LIPIŃSKIEGO (RULE OF 5)")
print("="*70)
# Kryteria Lipińskiego:
# 1. Masa molowa ≤ 500 Da
# 2. LogP ≤ 5
# 3. HBD ≤ 5
# 4. HBA ≤ 10
df['Lipinski_MW'] = df['MolWt'] <= 500
df['Lipinski_LogP'] = df['LogP'] <= 5
df['Lipinski_HBD'] = df['HBD'] <= 5
df['Lipinski_HBA'] = df['HBA'] <= 10
df['Lipinski_violations'] = (~df['Lipinski_MW']).astype(int) + \
(~df['Lipinski_LogP']).astype(int) + \
(~df['Lipinski_HBD']).astype(int) + \
(~df['Lipinski_HBA']).astype(int)
df['Lipinski_compliant'] = df['Lipinski_violations'] <= 1 # 0 lub 1 naruszenie = OK
# Statystyki
compliant = df['Lipinski_compliant'].sum()
non_compliant = (~df['Lipinski_compliant']).sum()
print(f"\n📊 Zgodność z regułą Lipińskiego:")
print(f" • Zgodne (≤1 naruszenie): {compliant:,} ({compliant/len(df)*100:.1f}%)")
print(f" • Niezgodne (>1 naruszenie): {non_compliant:,} ({non_compliant/len(df)*100:.1f}%)")
print(f"\n📈 Rozkład naruszeń:")
violations_dist = df['Lipinski_violations'].value_counts().sort_index()
for count, freq in violations_dist.items():
pct = (freq / len(df)) * 100
bar = '█' * int(pct / 2)
print(f" {count:.0f} naruszeń: {freq:5,} ({pct:5.2f}%) {bar}")
print(f"\n📋 Szczegółowe naruszenia:")
print(f" • MW > 500: {(~df['Lipinski_MW']).sum():,}")
print(f" • LogP > 5: {(~df['Lipinski_LogP']).sum():,}")
print(f" • HBD > 5: {(~df['Lipinski_HBD']).sum():,}")
print(f" • HBA > 10: {(~df['Lipinski_HBA']).sum():,}")
# Wizualizacja
fig, axes = plt.subplots(1, 3, figsize=(18, 5))
# 1. Pie chart - compliance
axes[0].pie([compliant, non_compliant],
labels=['Zgodne\n(Ro5)', 'Niezgodne'],
autopct='%1.1f%%', colors=['#2ecc71', '#e74c3c'],
startangle=90, textprops={'fontsize': 12, 'fontweight': 'bold'})
axes[0].set_title('Zgodność z regułą Lipińskiego', fontsize=13, fontweight='bold')
# 2. Barplot naruszeń
axes[1].bar(violations_dist.index, violations_dist.values, color='#9b59b6',
edgecolor='black', alpha=0.8, linewidth=1)
axes[1].set_xlabel('Liczba naruszeń Ro5', fontsize=12, fontweight='bold')
axes[1].set_ylabel('Liczba molekuł', fontsize=12, fontweight='bold')
axes[1].set_title('Rozkład naruszeń Ro5', fontsize=13, fontweight='bold')
axes[1].grid(axis='y', alpha=0.3)
# 3. Szczegółowe naruszenia
violation_types = ['MW > 500', 'LogP > 5', 'HBD > 5', 'HBA > 10']
violation_counts = [
(~df['Lipinski_MW']).sum(),
(~df['Lipinski_LogP']).sum(),
(~df['Lipinski_HBD']).sum(),
(~df['Lipinski_HBA']).sum()
]
axes[2].barh(range(len(violation_types)), violation_counts, color='#e74c3c',
edgecolor='black', alpha=0.8, linewidth=1)
axes[2].set_yticks(range(len(violation_types)))
axes[2].set_yticklabels(violation_types)
axes[2].set_xlabel('Liczba molekuł', fontsize=12, fontweight='bold')
axes[2].set_title('Typy naruszeń Ro5', fontsize=13, fontweight='bold')
axes[2].grid(axis='x', alpha=0.3)
for i, val in enumerate(violation_counts):
axes[2].text(val + 5, i, f'{val}', va='center', fontweight='bold')
plt.tight_layout()
plt.savefig('qm9_lipinski_ro5.png', dpi=300, bbox_inches='tight')
print("\n✅ Zapisano: qm9_lipinski_ro5.png")
plt.show()
# ============================================================
# 6. KORELACJE MIĘDZY DESKRYPTORAMI
# ============================================================
print("\n" + "="*70)
print("🔗 KORELACJE MIĘDZY DESKRYPTORAMI")
print("="*70)
# Wybór deskryptorów do analizy
descriptors_to_correlate = ['MolWt', 'LogP', 'TPSA', 'HBD', 'HBA', 'MolMR',
'HeavyAtoms', 'rotatable_bonds', 'frac_csp3']
df_corr = df[descriptors_to_correlate].copy()
# Macierz korelacji
corr_matrix = df_corr.corr()
print("\n📊 Najsilniejsze korelacje:")
# Ekstakcja korelacji (bez diagonali)
correlations_list = []
for i in range(len(corr_matrix.columns)):
for j in range(i+1, len(corr_matrix.columns)):
correlations_list.append({
'Prop1': corr_matrix.columns[i],
'Prop2': corr_matrix.columns[j],
'Correlation': corr_matrix.iloc[i, j]
})
df_corr_list = pd.DataFrame(correlations_list)
df_corr_list = df_corr_list.sort_values('Correlation', key=abs, ascending=False)
print("\nTop 10 najsilniejszych korelacji:")
for idx, row in df_corr_list.head(10).iterrows():
print(f" {row['Prop1']:20s} ↔ {row['Prop2']:20s} r = {row['Correlation']:+.3f}")
# Wizualizacja
fig, axes = plt.subplots(1, 2, figsize=(16, 6))
# 1. Heatmap
import seaborn as sns
sns.heatmap(corr_matrix, annot=True, fmt='.2f', cmap='RdBu_r',
center=0, square=True, linewidths=0.5, cbar_kws={"shrink": 0.8},
vmin=-1, vmax=1, ax=axes[0])
axes[0].set_title('Macierz korelacji deskryptorów', fontsize=13, fontweight='bold')
# 2. Top korelacje - barplot
top10 = df_corr_list.head(10).copy()
top10['Pair'] = top10['Prop1'] + ' ↔\n' + top10['Prop2']
colors_corr = ['#d73027' if x < 0 else '#4575b4' for x in top10['Correlation']]
axes[1].barh(range(len(top10)), top10['Correlation'], color=colors_corr,
edgecolor='black', alpha=0.8, linewidth=0.7)
axes[1].set_yticks(range(len(top10)))
axes[1].set_yticklabels(top10['Pair'], fontsize=8)
axes[1].set_xlabel('Współczynnik korelacji', fontsize=12, fontweight='bold')
axes[1].set_title('Top 10 najsilniejszych korelacji', fontsize=13, fontweight='bold')
axes[1].axvline(0, color='black', linewidth=0.8)
axes[1].grid(axis='x', alpha=0.3)
axes[1].invert_yaxis()
for i, val in enumerate(top10['Correlation']):
axes[1].text(val + 0.02 if val > 0 else val - 0.02, i,
f'{val:.2f}', va='center', ha='left' if val > 0 else 'right',
fontweight='bold', fontsize=8)
plt.tight_layout()
plt.savefig('qm9_descriptors_correlation.png', dpi=300, bbox_inches='tight')
print("\n✅ Zapisano: qm9_descriptors_correlation.png")
plt.show()
# ============================================================
# 7. PODSUMOWANIE
# ============================================================
print("\n" + "="*70)
print("📋 PODSUMOWANIE WŁAŚCIWOŚCI FIZYKOCHEMICZNYCH")
print("="*70)
print(f"""
💧 LIPOFILOWOŚĆ (LogP):
• Średnia: {df['LogP'].mean():.3f}
• Zakres: {df['LogP'].min():.2f} - {df['LogP'].max():.2f}
• Hydrofilowe (<1): {(df['LogP'] < 1).sum():,} ({(df['LogP'] < 1).sum()/len(df)*100:.1f}%)
• Lipofilowe (≥3): {(df['LogP'] >= 3).sum():,} ({(df['LogP'] >= 3).sum()/len(df)*100:.1f}%)
🧊 TPSA (Polar Surface Area):
• Średnia: {df['TPSA'].mean():.2f} Ų
• BBB permeable: {bbb_permeable:,} ({bbb_permeable/len(df)*100:.1f}%)
• CNS optimal: {cns_optimal:,} ({cns_optimal/len(df)*100:.1f}%)
🔗 WIĄZANIA WODOROWE:
• Średnia HBD: {df['HBD'].mean():.2f}
• Średnia HBA: {df['HBA'].mean():.2f}
• Max HBD: {df['HBD'].max():.0f}
• Max HBA: {df['HBA'].max():.0f}
💊 REGUŁA LIPIŃSKIEGO (Ro5):
• Zgodne molekuły: {compliant:,} ({compliant/len(df)*100:.1f}%)
• 0 naruszeń: {(df['Lipinski_violations']==0).sum():,}
• 1 naruszenie: {(df['Lipinski_violations']==1).sum():,}
• >1 naruszeń: {(df['Lipinski_violations']>1).sum():,}
🔗 KORELACJE:
• Najsilniejsza: {df_corr_list.iloc[0]['Prop1']} ↔ {df_corr_list.iloc[0]['Prop2']} (r={df_corr_list.iloc[0]['Correlation']:.3f})
• Silne korelacje (|r|>0.7): {(df_corr_list['Correlation'].abs() > 0.7).sum()}
✅ WNIOSKI:
• QM9 zawiera głównie małe, relatywnie hydrofilowe molekuły
• Wysoka zgodność z regułą Lipińskiego ({compliant/len(df)*100:.0f}%)
• Niska TPSA → dobra przepuszczalność błon
• Deskryptory silnie skorelowane z masą molową
""")
print("="*70)
print("✅ ANALIZA WŁAŚCIWOŚCI FIZYKOCHEMICZNYCH ZAKOŃCZONA!")
print("="*70)
======================================================================
WŁAŚCIWOŚCI FIZYKOCHEMICZNE - QM9 DATASET
======================================================================
⚗️ ETAP 1: OBLICZANIE DESKRYPTORÓW FIZYKOCHEMICZNYCH
----------------------------------------------------------------------
Obliczam właściwości fizykochemiczne...
✅ Obliczono 7 właściwości dla 1000 molekuł
📊 Statystyki właściwości fizykochemicznych:
LogP MolMR TPSA HBD HBA \
count 1000.000000 1000.000000 1000.000000 1000.000000 1000.000000
mean 0.104381 31.821830 37.802210 0.972000 2.053000
std 1.173836 3.606108 21.399043 0.889948 1.093337
min -3.467400 14.868400 0.000000 0.000000 0.000000
25% -0.660365 29.675525 21.510000 0.000000 1.000000
50% 0.158600 31.994000 36.510000 1.000000 2.000000
75% 0.923925 34.185250 52.512500 2.000000 3.000000
max 3.222700 41.365000 112.360000 4.000000 5.000000
LabuteASA BalabanJ
count 1000.000000 1000.000000
mean 52.579936 2.448441
std 3.722869 0.525115
min 24.605606 0.000000
25% 51.782555 2.071725
50% 53.587200 2.269522
75% 54.820772 2.693324
max 58.088080 4.837063
======================================================================
💧 ANALIZA LogP (LIPOFILOWOŚĆ)
======================================================================
📊 Statystyki LogP:
• Średnia: 0.104
• Mediana: 0.159
• Std: 1.174
• Min: -3.467
• Max: 3.223
📈 Percentyle LogP:
• 10%: -1.387
• 25%: -0.660
• 50%: 0.159
• 75%: 0.924
• 90%: 1.561
📋 Klasyfikacja lipofilowości:
• Bardzo hydrofilowe (LogP < -1): 179 (17.9%)
• Hydrofilowe (-1 ≤ LogP < 1): 592 (59.2%)
• Umiarkowane (1 ≤ LogP < 3): 227 (22.7%)
• Lipofilowe (3 ≤ LogP < 5): 2 (0.2%)
• Bardzo lipofilowe (LogP ≥ 5): 0 (0.0%)
⚠️ Outliers LogP (metoda IQR):
• Liczba outliers: 7 (0.70%)
• Zakres normalny: [-3.04, 3.30]
C:\Users\slast\AppData\Local\Temp\ipykernel_6436\1039049785.py:114: MatplotlibDeprecationWarning: The 'labels' parameter of boxplot() has been renamed 'tick_labels' since Matplotlib 3.9; support for the old name will be dropped in 3.11. bp = axes[0, 1].boxplot([df['LogP'].dropna()], labels=['LogP'],
✅ Zapisano: qm9_logp_analysis.png
====================================================================== 🧊 ANALIZA TPSA (TOPOLOGICAL POLAR SURFACE AREA) ====================================================================== 📊 Statystyki TPSA: • Średnia: 37.80 Ų • Mediana: 36.51 Ų • Std: 21.40 Ų • Max: 112.36 Ų 💊 Drug-likeness (TPSA): • BBB permeable (< 140 Ų): 1,000 (100.0%) • CNS optimal (< 90 Ų): 991 (99.1%) ✅ Zapisano: qm9_tpsa_analysis.png
====================================================================== 🔗 ANALIZA WIĄZAŃ WODOROWYCH (HBD/HBA) ====================================================================== 📊 Statystyki: • Średnia HBD (donory): 0.97 • Średnia HBA (akceptory): 2.05 • Max HBD: 4 • Max HBA: 5 📈 Rozkład HBD: 0 donorów: 348 (34.80%) ███████████ 1 donorów: 392 (39.20%) █████████████ 2 donorów: 202 (20.20%) ██████ 3 donorów: 56 ( 5.60%) █ 📈 Rozkład HBA: 0 akceptorów: 66 ( 6.60%) ██ 1 akceptorów: 249 (24.90%) ████████ 2 akceptorów: 355 (35.50%) ███████████ 3 akceptorów: 247 (24.70%) ████████ 4 akceptorów: 62 ( 6.20%) ██ 5 akceptorów: 21 ( 2.10%) ✅ Zapisano: qm9_hbond_analysis.png
====================================================================== 💊 REGUŁA LIPIŃSKIEGO (RULE OF 5) ====================================================================== 📊 Zgodność z regułą Lipińskiego: • Zgodne (≤1 naruszenie): 1,000 (100.0%) • Niezgodne (>1 naruszenie): 0 (0.0%) 📈 Rozkład naruszeń: 0 naruszeń: 1,000 (100.00%) ██████████████████████████████████████████████████ 📋 Szczegółowe naruszenia: • MW > 500: 0 • LogP > 5: 0 • HBD > 5: 0 • HBA > 10: 0
C:\Users\slast\AppData\Local\Temp\ipykernel_6436\1039049785.py:358: UserWarning: Tight layout not applied. The left and right margins cannot be made large enough to accommodate all Axes decorations. plt.tight_layout()
✅ Zapisano: qm9_lipinski_ro5.png
====================================================================== 🔗 KORELACJE MIĘDZY DESKRYPTORAMI ====================================================================== 📊 Najsilniejsze korelacje: Top 10 najsilniejszych korelacji: MolWt ↔ HeavyAtoms r = +0.927 TPSA ↔ HBA r = +0.821 LogP ↔ HBD r = -0.638 MolMR ↔ HeavyAtoms r = +0.578 TPSA ↔ HBD r = +0.572 LogP ↔ MolMR r = +0.495 LogP ↔ TPSA r = -0.476 MolWt ↔ MolMR r = +0.470 HBA ↔ MolMR r = -0.462 TPSA ↔ frac_csp3 r = -0.447 ✅ Zapisano: qm9_descriptors_correlation.png
====================================================================== 📋 PODSUMOWANIE WŁAŚCIWOŚCI FIZYKOCHEMICZNYCH ====================================================================== 💧 LIPOFILOWOŚĆ (LogP): • Średnia: 0.104 • Zakres: -3.47 - 3.22 • Hydrofilowe (<1): 771 (77.1%) • Lipofilowe (≥3): 2 (0.2%) 🧊 TPSA (Polar Surface Area): • Średnia: 37.80 Ų • BBB permeable: 1,000 (100.0%) • CNS optimal: 991 (99.1%) 🔗 WIĄZANIA WODOROWE: • Średnia HBD: 0.97 • Średnia HBA: 2.05 • Max HBD: 4 • Max HBA: 5 💊 REGUŁA LIPIŃSKIEGO (Ro5): • Zgodne molekuły: 1,000 (100.0%) • 0 naruszeń: 1,000 • 1 naruszenie: 0 • >1 naruszeń: 0 🔗 KORELACJE: • Najsilniejsza: MolWt ↔ HeavyAtoms (r=0.927) • Silne korelacje (|r|>0.7): 2 ✅ WNIOSKI: • QM9 zawiera głównie małe, relatywnie hydrofilowe molekuły • Wysoka zgodność z regułą Lipińskiego (100%) • Niska TPSA → dobra przepuszczalność błon • Deskryptory silnie skorelowane z masą molową ====================================================================== ✅ ANALIZA WŁAŚCIWOŚCI FIZYKOCHEMICZNYCH ZAKOŃCZONA! ======================================================================